pythonStagYY/P1_figS5.py at main · chingchenn/pythonStagYY · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 15 18:55:19 2024

@author: chingchen
"""

import pandas as pd
import numpy as np
import numpy.ma as ma
# from matplotlib import cm
# import matplotlib  as mpl
# from scipy.misc import derivative
import matplotlib.pyplot as plt

labelsize = 20
bwith = 3

### PATH ###
path = '/Users/chingchen/Desktop/data/'
workpath = '/Users/chingchen/Desktop/StagYY_Works/thermal_evolution_v2/'
modelpath = '/Users/chingchen/Desktop/model/'
figpath = '/Users/chingchen/Desktop/figure/'
colors=['#282130','#3CB371','#4682B4','#CD5C5C','#97795D','#414F67','#4198B9','#3CB371']
header_list = ['time_Gyr','Prad','Ptidal','Fcore','Pint','Hint','conv',
               'melt','P','zbot','%vol','Tbot','Tm','Fbot','Ftop','dlid','T_core']


fig1,(aa1,aa2) = plt.subplots(2,2,figsize=(17,11))
ax1 = aa1[0]
ax2 = aa1[1]
ax3 = aa2[0]
ax4 = aa2[1]

label_list=['0.6 TW','0.8 TW','1.0 TW','1.2 TW']
# ---------------------------------------- figure 3 dash 1.0 TW --------------------------------
model_list = ['Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P0.6TW_1.5wt%_D5.0km-NH3', # power
              'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P0.8TW_1.5wt%_D5.0km-NH3',
              'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
              'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.2TW_1.5wt%_D5.0km-NH3',]

for i, model in enumerate(model_list):
    data = pd.read_csv(workpath+model+'_Hvar_2_thermal-evolution.dat',
        header=None,names=header_list,delim_whitespace=True)[2:].reset_index(drop=True)
    data = data.replace('D','e',regex=True).astype(float) # data convert to float
    x = data.time_Gyr
    mask_cond = data.conv
    mask_conv = ~ma.array(data.Ftop, mask = data.conv).mask
    zbot_cond = ma.array(data.Ftop, mask = mask_cond)
    zbot_conv = ma.array(data.Ftop, mask = mask_conv)
    ax1.plot(x,zbot_conv,color=colors[i],lw=3,label=label_list[i])
    ax1.plot(x,zbot_cond,color='orange')


label_list=['10$^{13}$','10$^{13.5}$','10$^{13.75}$','10$^{14}$']
# ---------------------------------------- figure 3 dash 1.0 TW --------------------------------
model_list = ['Europa-tidal5_period0.14Gyr_emx10%_eta1.0d13_P1.0TW_1.5wt%_D5.0km-NH3',# viscosity
             'Europa-tidal5_period0.14Gyr_emx10%_eta3.2d13_P1.0TW_1.5wt%_D5.0km-NH3',
             'Europa-tidal5_period0.14Gyr_emx10%_eta5.6d13_P1.0TW_1.5wt%_D5.0km-NH3',
             'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
             # 'Europa-tidal5_period0.14Gyr_emx10%_eta3.2d14_P1.0TW_1.5wt%_D5.0km-NH3',
             ]
for i, model in enumerate(model_list):
    data = pd.read_csv(workpath+model+'_Hvar_2_thermal-evolution.dat',
        header=None,names=header_list,delim_whitespace=True)[2:].reset_index(drop=True)
    data = data.replace('D','e',regex=True).astype(float) # data convert to float
    x = data.time_Gyr
    mask_cond = data.conv
    mask_conv = ~ma.array(data.Ftop, mask = data.conv).mask
    zbot_cond = ma.array(data.Ftop, mask = mask_cond)
    zbot_conv = ma.array(data.Ftop, mask = mask_conv)
    ax2.plot(x,zbot_conv,color=colors[i],lw=3,label=label_list[i])
    ax2.plot(x,zbot_cond,color='orange')


label_list=['vol = 0.%','vol = 1.5%','vol = 3.0%']
# ---------------------------------------- figure 4 dashed 1.0 TW--------------------------------
model_list = ['Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.0TW_0.0wt%_D5.0km-NH3',# composition
              'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
              'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.0TW_3.0wt%_D5.0km-NH3',
              ]
for i, model in enumerate(model_list):
    data = pd.read_csv(workpath+model+'_Hvar_2_thermal-evolution.dat',
        header=None,names=header_list,delim_whitespace=True)[2:].reset_index(drop=True)
    data = data.replace('D','e',regex=True).astype(float) # data convert to float
    x = data.time_Gyr
    mask_cond = data.conv
    mask_conv = ~ma.array(data.Ftop, mask = data.conv).mask
    zbot_cond = ma.array(data.Ftop, mask = mask_cond)
    zbot_conv = ma.array(data.Ftop, mask = mask_conv)
    ax3.plot(x,zbot_conv,color=colors[i],lw=2,label=label_list[i])
    ax3.plot(x,zbot_cond,color='orange')
label_list = ['5%','20%','50%']
# ---------------------------------------- figure 3 dash 1.0 TW --------------------------------
model_list = ['Europa-tidal5_period0.14Gyr_emx5%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
              # 'Europa-tidal5_period0.14Gyr_emx10%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
              'Europa-tidal5_period0.14Gyr_emx20%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
              'Europa-tidal5_period0.14Gyr_emx50%_eta1.0d14_P1.0TW_1.5wt%_D5.0km-NH3',
              ]
for i, model in enumerate(model_list):
    data = pd.read_csv(workpath+model+'_Hvar_2_thermal-evolution.dat',
        header=None,names=header_list,delim_whitespace=True)[2:].reset_index(drop=True)
    data = data.replace('D','e',regex=True).astype(float) # data convert to float
    x = data.time_Gyr
    mask_cond = data.conv
    mask_conv = ~ma.array(data.Ftop, mask = data.conv).mask
    zbot_cond = ma.array(data.Ftop, mask = mask_cond)
    zbot_conv = ma.array(data.Ftop, mask = mask_conv)
    ax4.plot(x,zbot_conv,color=colors[i],lw=3,label=label_list[i])
    ax4.plot(x,zbot_cond,color='orange')
#  ------------------------------ figure setting ------------------------------
ax3.set_xlabel('Time (Gyr)',fontsize=labelsize)
ax4.set_xlabel('Time (Gyr)',fontsize=labelsize)
for aa in [ax1,ax2,ax3,ax4]:
    aa.legend(fontsize=labelsize)
    aa.set_ylim(0,100)
    aa.set_ylabel('Surface heat flux',fontsize = labelsize)
    aa.minorticks_on()
    aa.tick_params(which='minor', length=5, width=2, direction='in')
    aa.tick_params(labelsize=labelsize,width=3,length=10,right=True, top=True,direction='in',pad=15)
    aa.set_xlim(0,4.55)
    aa.grid()
    for axis in ['top','bottom','left','right']:
        aa.spines[axis].set_linewidth(bwith)
# fig1.savefig('/Users/chingchen/Desktop/StagYY_Works/paper_europa_ice_shell/figureS5_v7.pdf')