-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsfr.py
233 lines (198 loc) · 10.7 KB
/
sfr.py
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
"""
Loads star formation rate data from the file with format:
{"z", "tStart", "tEnd", "totStarMass", "totPop3StarMass",
"totPollStarMass", "totPrimordStarMass", "totGasMass",
"totPristGasMass", "totSubcritStarMass", "totNonPrimordStarMass"}
Data is accessed by sfrd, sfrdP3, sfrdCP3 fields
"""
import numpy as np
from scipy.interpolate import interp1d
from astropy.cosmology import FlatLambdaCDM
from astropy import units as u
import astropy
class SFR:
"""
Uses data generated by the extract_star_data class to compute the
star formation rate density (SFRD). User specifies the simulation
boxsize (comoving Mpc/h), the path to the star particle data and
the value for h = H0/100.
size : float
Boxsize in Mpc/h
prefix : str ("./")
path to the star particle data summary file
hubble : float (0.71)
H0/(100 km/s/Mpc)
"""
def __init__(self, size, prefix="./", hubble=0.71):
self.prefix = prefix
self.hub = hubble
self.size = size
self.cosmo = FlatLambdaCDM(H0=hubble * 100.0,
Om0=0.267,
Ob0=0.0449,
name="myCosmo")
self.cm_vol = (self.size / self.hub)**3
return
def load(self, file="StarParticle-data3Mpc.txt"):
"""
Loads data expecting the following columns:
z tStart tEnd totStarMass totPop3StarMass totPollStarMass totPrimordStarMass totGasMass \
totPristGasMass totSubcritStarMass totNonPrimordStarMass
User specified the filename to read.
Once the data is loaded, the user can access the fields below:
redshift : an array of redshifts at which the SFRD value is computed
sfrdP3 : the SFRD for Pop III stars
sfrdCP3 : the SFRD for Pop III stars that does not consider subgrid
mixing
"""
self.file = file
self.rawSFRdata = np.loadtxt(self.prefix + self.file, skiprows=1)
# Determine delta_t between each entry (each entry is for a single time/redshift)
# Use the latest time as the time associated with the SFRD up to that time
# Compute delta_mass
# Compute average SFRD over the preceeding interval: delta_mass/delta_t * nmlze
# where nmlze is the comoving volume of the simulation
self.time = self.rawSFRdata[:, 2] # Get all end times
self.delta_t = np.diff(self.time) # Time between entries
# End time point. Still calling it "mid_time"
self.mid_time = self.time[1:]
# How much solar mass was created over the previous interval?
self.delta_m = np.diff(self.rawSFRdata[:, 3])
# If we have little or no star formation, but SN, dM could be negative!!
self.delta_m[self.delta_m < 0.0] = np.nan
self.sfrd = self.delta_m / self.delta_t / self.cm_vol
# self.sfrd[np.isnan(self.sfrd)] = 1e-10 # No star formation during delta t, set to small #
# Compute the redshift associated with the time...
self.midz = []
for t in self.mid_time:
# Using astropy is ok, although it doesn't precisely match the orig file's redshifts
self.midz.append(
astropy.cosmology.z_at_value(self.cosmo.age, t * u.yr))
# Associate the delta sfrd with the later z
self.redshift = self.rawSFRdata[1:, 0]
self.delta_mP3 = np.diff(self.rawSFRdata[:, 4])
self.delta_mP3[self.delta_mP3 < 0.0] = np.nan
self.sfrdP3 = self.delta_mP3 / self.delta_t / self.cm_vol
# self.sfrdP3[np.isnan(self.sfrdP3)] = 1e-10
self.delta_mCP3 = np.diff(self.rawSFRdata[:, 9])
self.delta_mCP3[self.delta_mCP3 < 0.0] = np.nan
self.sfrdCP3 = self.delta_mCP3 / self.delta_t / self.cm_vol
# self.sfrdCP3[np.isnan(self.sfrdCP3)] = 1e-10
# Experimental ... use fitting to interpolate the values and compute the
# SFRD at arbitrary redshifts...
# def set_baseline(self, base):
# self.baseline={}
## self.baseline['tot'] = interp1d(base.zData,base.totStarMass,kind='linear')
## self.baseline['pop3'] = interp1d(base.zData,base.popIIIMass, kind='linear')
## self.baseline['cpop3'] = interp1d(base.zData,base.subCritMass,kind='linear')
# def get_diff_funcs(self,newVal,base=None):
# if not 'self.baseline' in locals and base==None:
## print("Specify baseline for differencing... ")
# return
# if not 'self.baseline' in locals:
# self.initBaseline(base)
# self.other={}
## self.other['tot'] = interp1d(newVal.zData,newVal.totStarMass,kind='linear')
## self.other['pop3'] = interp1d(newVal.zData,newVal.popIIIMass, kind='linear')
## self.other['cpop3'] = interp1d(newVal.zData,newVal.subCritMass,kind='linear')
## zmax = min(newVal.zData.max(),self.baseline.zData.max())
## zmin = max(newVal.zData.min(),self.baseline.zData.min())
## theRange = np.linspace(zmin,zmax,10)
# return theRange,interp1d(theRange, -self.baseline['tot'](theRange) + self.other['tot'](theRange),kind='linear'), \
## interp1d(theRange, -self.baseline['pop3'](theRange) + self.other['pop3'](theRange),kind='linear'), \
## interp1d(theRange, -self.baseline['cpop3'](theRange) + self.other['cpop3'](theRange), kind='linear')
def loadNew(self, file):
"""
Loads data expecting the following columns:
z time tStart tEnd totStarMass totPop3StarMass totPollStarMass totPrimordStarMass totGasMass \
totPristGasMass totSubcritStarMass totNonPrimordStarMass
User specified the filename to read.
Once the data is loaded, the user can access the fields below:
redshift : an array of redshifts at which the SFRD value is computed
sfrdP3 : the SFRD for Pop III stars
sfrdCP3 : the SFRD for Pop III stars that does not consider subgrid
mixing
"""
self.file = file
self.rawSFRdata = np.loadtxt(self.prefix + self.file, skiprows=1)
# Determine delta_t between each entry (each entry is for a single time/redshift)
# Use the time col as the time associated with the SFRD up to that time
# Compute delta_mass
# Compute average SFRD over the preceeding interval: delta_mass/delta_t * nmlze
# where nmlze is the comoving volume of the simulation
self.time = self.rawSFRdata[:, 1] # Get all end times
self.delta_t = np.diff(self.time) # Time between entries
# End time point. Still calling it "mid_time"
self.mid_time = self.time[1:]
# How much solar mass was created over the previous interval?
self.delta_m = np.diff(self.rawSFRdata[:, 4])
# If we have little or no star formation, but SN, dM could be negative!!
self.delta_m[self.delta_m < 0.0] = np.nan
self.sfrd = self.delta_m / (self.delta_t * 1e6) / self.cm_vol
# self.sfrd[np.isnan(self.sfrd)] = 1e-10 # No star formation during delta t, set to small #
# Compute the redshift associated with the time...
self.midz = []
for t in self.mid_time:
# Using astropy is ok, although it doesn't precisely match the orig file's redshifts
self.midz.append(
astropy.cosmology.z_at_value(self.cosmo.age,
t * u.Myr)) # New time is Myr
# Associate the delta sfrd with the later z
self.redshift = self.rawSFRdata[1:, 0]
self.delta_mP3 = np.diff(self.rawSFRdata[:, 5])
self.delta_mP3[self.delta_mP3 < 0.0] = np.nan
self.sfrdP3 = self.delta_mP3 / (self.delta_t * 1e6) / self.cm_vol
# self.sfrdP3[np.isnan(self.sfrdP3)] = 1e-10
self.delta_mCP3 = np.diff(self.rawSFRdata[:, 10])
self.delta_mCP3[self.delta_mCP3 < 0.0] = np.nan
self.sfrdCP3 = self.delta_mCP3 / (self.delta_t * 1e6) / self.cm_vol
# self.sfrdCP3[np.isnan(self.sfrdCP3)] = 1e-10
def loadnf(self, file="starData.txt"):
"""
Loads data expecting the following columns:
z time toSPM totPIIIM totPIIM totCPIIIM totZPM
User specifies the filename to read.
file : (starData.txt) the file to load
Once the data is loaded, the user can access the fields below:
redshift : an array of redshifts at which the SFRD value is computed
time : the proper time associated with z, accounting for
my standard cosmology
sfrdP3 : the SFRD for Pop III stars
sfrdCP3 : the SFRD for Pop III stars that does not consider subgrid
mixing
"""
self.file = file
self.sfrFile = np.loadtxt(self.prefix + self.file, skiprows=1)
# Determine delta_t between each entry (each entry is for a single time/redshift)
# Use the latest time as the time associated with the SFRD up to that time
# Compute delta_mass
# Compute average SFRD over the preceeding interval: delta_mass/delta_t * nmlze
# where nmlze is the comoving volume of the simulation
# z time totSPM totPIIIM totPIIM totCPIIIM totZPM
# 0 1 2 3 4 5 6
self.time = self.sfrFile[:, 1] # Get time in Myr
self.delta_t = np.diff(self.time) * 1e6 # Time between entries, yr
# End time point. Still calling it "mid_time"
# Ignore the first dt since it isn't correct
self.mid_time = self.time[1:]
# How much solar mass was created over the previous interval?
self.delta_m = np.diff(self.sfrFile[:, 2])
# If we have little or no star formation, but SN, dM could be negative!!
self.delta_m[self.delta_m < 0.0] = np.nan
self.sfrd = self.delta_m / self.delta_t / self.cm_vol
# self.sfrd[np.isnan(self.sfrd)] = 1e-10 # No star formation during delta t, set to small #
# Compute the redshift associated with the time...
self.midz = []
for t in self.mid_time:
# Using astropy is ok, although it doesn't precisely match the orig file's redshifts
self.midz.append(
astropy.cosmology.z_at_value(self.cosmo.age, t * u.Myr))
# Associate the delta sfrd with the later z
self.redshift = self.sfrFile[1:, 0]
self.delta_mP3 = np.diff(self.sfrFile[:, 3])
self.delta_mP3[self.delta_mP3 < 0.0] = np.nan
self.sfrdP3 = self.delta_mP3 / self.delta_t / self.cm_vol
# self.sfrdP3[np.isnan(self.sfrdP3)] = 1e-10
self.delta_mCP3 = np.diff(self.sfrFile[:, 5])
self.delta_mCP3[self.delta_mCP3 < 0.0] = np.nan
self.sfrdCP3 = self.delta_mCP3 / self.delta_t / self.cm_vol