-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathdatalooper.py
More file actions
174 lines (157 loc) · 6.01 KB
/
datalooper.py
File metadata and controls
174 lines (157 loc) · 6.01 KB
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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
'''
Created 12.08.2022, author @piiamt
LOOPS THROUGH RAN SIMULATIONS AND RUN IDLS AND THEN OVERVIEW
PLOTTING, SAVING THE OVERVIEW PLOT AND NPY FILES
'''
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sys
import os, os.path
from os.path import expanduser
home = '/projects/astro3/nobackup/piia/'
import shutil
from idlpy import *
plt.rcParams.update({'font.size' : 10,
'mathtext.fontset' : 'cm',
'font.family' : 'serif',
'xtick.direction' :'in',
'ytick.direction' :'in',
#'axes.grid' :True,
'grid.alpha' : 0.7,
'grid.linestyle' : 'dashed'
})
colors = np.array(['k', '#009e73', '#d55e00', '#cc79a7', '#0072b2', '#e69f00', '#56b4e9'])
RGBc = np.array([[0,0,0,0.2],[0,158/255,112/255,0.2],[213/255,94/255,0.0,0.2],[204/255,121/255,167/255,0.2],
[0,114/255,178/255,0.2],[230/255,159/255,0,0.2],[86/255,180/255,233/255,0.2]])
# aka black green vermillion pink blue orange lskdfs
##### CONSTANTS #####
apla_min = 0.7
apla_max = 0.85
apla_inc = 0.05
mpla_min = 0.9
mpla_max = 1.05
mpla_inc = 0.1
myr = 365.25*24*60*60*1e6 # seconds
mearth = 5.9722e24 # kg
msun = 1.98847e30 # kg
AU = 1.496e11 # m
R0 = 250*1000 # meters
rho0 = 1.421573328e20 / ((4/3)*np.pi*R0**3) # kg/m^3
M0 = rho0 * (4/3)*np.pi*R0**3 # kg (M of initial planetesimal)
Aplas = np.arange(apla_min, apla_max, apla_inc)
Mplas = np.arange(mpla_min, mpla_max, mpla_inc)
def runidls(sim):
'''
With IDL pre-opened, runs adap_read_var, obj=var, /lsave,
saving only the end results of the simulation
to preserve memory and then adap_read_ts, obj=ts, /lsave
'''
os.environ['SHELL'] = 'tcsh'
os.chdir(home+sim)
# open IDL
# do "adap_read_var, /lsave"
# do "adap_read_ts, /lsave"
# close IDL
def runoverview(sim):
os.chdir(home+sim)
timefile = home+sim+'/time_series.idl'
timef = readsav(timefile)
tser = list(timef.values())[0][0]
### plotting filter:
t5 = np.where(tser.t>(5*Myr_to_s))
ts = tser.t[t5]/Myr_to_s
#make the plot
fig, ax = plt.subplots(3,4, figsize=(12,9), sharex=True,
gridspec_kw={'wspace':0.2,'hspace':0.2})
# plt.gca().set_yscale('log')
fig.suptitle(folder, fontsize=15)
ax[0,0].set_ylabel(r'Mass ($M_{\oplus}$)')
ax[0,0].set_yscale('log')
ax[0,0].plot(ts, tser.matm[t5]/mearth, c=colors[1])
ax[0,0].set_title('Atmosphere M', fontsize=10)
ax[0,1].plot(ts, tser.mpla[t5]/mearth, c=colors[3])
ax[0,1].set_ylim(np.mean(tser.mpla[t5])/mearth-0.005,
np.mean(tser.mpla[t5]/mearth+0.005))
ax[0,1].yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[0,1].set_title('Planet M', fontsize=10)
ax[0,2].plot(ts, tser.matm_o2[t5]/tser.matm[t5], c=colors[0])
ax[0,2].set_ylim(-0.05,1.05)
ax[0,2].set_title(r'O$_2$ / Matm', fontsize=10)
ax[0,3].plot(ts, tser.matm_h2o[t5]/tser.matm[t5], c=colors[5])
ax[0,3].set_ylabel(r'Atm mass fraction')
ax[0,3].yaxis.set_label_position('right')
ax[0,3].set_ylim(-0.05,1.05)
ax[0,3].set_title(r'H$_2$O / Matm', fontsize=10)
ax[1,0].set_ylabel(r'Mass ($M_{\oplus}$)')
ax[1,0].set_yscale('log')
ax[1,0].plot(ts, tser.menv[t5]/mearth, c=colors[4])
ax[1,0].set_title('Envelope M', fontsize=10)
ax[1,1].plot(ts, tser.tsur[t5], c=colors[5])
ax[1,1].set_title('Surface T (K)', fontsize=10)
ax[1,2].plot(ts, tser.matm_co[t5]/tser.matm[t5], c=colors[2])
ax[1,2].set_ylim(-0.05,1.05)
ax[1,2].set_title(r'CO / Matm', fontsize=10)
ax[1,3].plot(ts, tser.matm_co2[t5]/tser.matm[t5], c=colors[1])
ax[1,3].set_ylabel(r'Atm mass fraction')
ax[1,3].yaxis.set_label_position('right')
ax[1,3].set_ylim(-0.05,1.05)
ax[1,3].set_title(r'CO$_2$ / Matm', fontsize=10)
ax[2,0].set_ylabel(r'Mass ($M_{\oplus}$)')
ax[2,0].set_yscale('log')
ax[2,0].plot(ts, tser.msolid[t5]/mearth, c=colors[0])
ax[2,0].set_title('Solids M', fontsize=10)
ax[2,1].plot(ts, tser.psur[t5]/1e5, c=colors[1])
ax[2,1].set_yscale('log')
ax[2,1].set_title('Surface pressure (bar)', fontsize=10)
ax[2,2].plot(ts, tser.matm_n2[t5]/tser.matm[t5], c=colors[3])
ax[2,2].set_ylim(-0.05,1.05)
ax[2,2].set_title(r'N$_2$ / Matm', fontsize=10)
ax[2,3].plot(ts, tser.matm_h2[t5]/tser.matm[t5], c=colors[4])
ax[2,3].set_ylim(-0.05,1.05)
ax[2,3].set_title(r'H$_2$ / Matm', fontsize=10)
ax[2,3].set_ylabel(r'Atm mass fraction')
ax[2,3].yaxis.set_label_position('right')
ax[2,0].set_xlabel('Time (Myr)')
ax[2,1].set_xlabel('Time (Myr)')
ax[2,2].set_xlabel('Time (Myr)')
ax[2,3].set_xlabel('Time (Myr)')
OGhome = expanduser('~')
plt.savefig(OGhome+'/python/autoplots'+folder+'.png',
bbox_inches='tight')
### Saving everything into NPY files for faster later plotting
### We are already in the simulation's dir so just saving:
dirname = home + folder + '/pysave'
os.mkdir(dirname) # The .npy files will be saved in this dir
np.save(dirname+'/t', tser.t)
np.save(dirname+'/mpla', tser.mpla)
np.save(dirname+'/menv', tser.menv)
np.save(dirname+'/matm', tser.matm)
np.save(dirname+'/msolid', tser.msolid)
np.save(dirname+'/psur', tser.psur)
np.save(dirname+'/tsur', tser.tsur)
np.save(dirname+'/matm_o2', tser.matm_o2)
np.save(dirname+'/matm_co', tser.matm_co)
np.save(dirname+'/matm_co2', tser.matm_co2)
np.save(dirname+'/matm_h2o', tser.matm_h2o)
np.save(dirname+'/matm_n2', tser.matm_n2)
np.save(dirname+'/matm_h2', tser.matm_h2)
### I will not save rpla bc that is in the variables.idl
### file and opening that is very lengthy
return(None)
### First we run the loop to IDL all the
### simulations that have been run
# OPEN IDL SOMEHOW
os.environ['SHELL'] = 'tcsh'
# SOMETHONG SOMETHING MAGIC
for apla in Aplas:
for mpla in Mplas:
sim = '/'+apla+'AU_'+mpla+'Mearth'
os.chdir(home+sim)
runidls(sim)
### Now we run the loop for overviewplotting
### all the data of the simulations
for apla in Aplas:
for mpla in Mplas:
sim = '/'+apla+'AU_'+mpla+'Mearth'
runoverview(sim)