-
Notifications
You must be signed in to change notification settings - Fork 0
/
buildFluxTables.py
238 lines (205 loc) · 10.9 KB
/
buildFluxTables.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
234
235
236
237
238
import numpy as np
from scipy.io import readsav
import re
#from astropy.cosmology import FlatLambdaCDM
from astropy import units as u
import astropy
from matplotlib.pylab import *
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
from matplotlib.ticker import NullFormatter
rcParams['figure.figsize'] = (10,8)
rcParams['font.size'] = 22
import sys
sys.path.append('/Users/earnric/Google Drive/ASU/Codes/PythonCode/modules')
import loadfilt as lf
import igm as lyA
import itertools
import os
import subprocess
import glob
import gc
import linecache
#
# Build the file from the flux data and redshift
#
def buildOutFile(z,allAges,ages,wavelns,LperA):
prevAge = 0.0
outfile = []
print("working {} at z={:.1f}".format(fname,z))
absorb = lyA.lyTauC(z) # Create a lyman forest absorption function for z... ***** NEW truncated absorption
for age in ages: # Once for each age group in the list
if age - prevAge < 0.025:
continue # We don't super-tight spacing in age...
print("log age={:.4f}".format(age))
prevAge = age
ageCond = (allAges == age) # select records based on log age field...
# Convert to freq & Lumin/s/Hz
freq = (wavelns[ageCond][::-1]).to(u.Hz, equivalencies=u.spectral()) # Reverse array order
LperHz = (LperA[ageCond] * wavelns[ageCond]**2/astropy.constants.c).to(u.erg/u.s/u.Hz)[::-1] # Reverse array order...
rsWaveln, rsLperA = lyA.rsSEDwavelen(wavelns[ageCond], LperA[ageCond], z) # Redshift SED in wavelen
rsFreq, rsLperHz = lyA.rsSEDfreq(freq, LperHz, z) # Redshift SED in freq
lyForFluxHz = (rsLperHz * absorb(rsWaveln[::-1]))
# Need to extract only the data for one sed (at a single age)
# to compute the flux in the filters
rsWavelnOneAge = rsWaveln
rsFreqOneAge = rsFreq
lyForFluxHzOneAge = lyForFluxHz
# Create first part [log age, z, ...]
# Flux is in erg/s/Hz/cm^2
jwstNormFlux = np.array([np.trapz(jwstFilters[aFilt](rsWavelnOneAge[::-1])/rsFreqOneAge,rsFreqOneAge).value for aFilt in jwstFilters])
jwstNormFlux[jwstNormFlux <= 0.0] = -1.0
jwstFlux = np.array([(np.trapz(jwstFilters[aFilt](rsWavelnOneAge[::-1])*lyForFluxHzOneAge/rsFreqOneAge,rsFreqOneAge)).value
for aFilt in jwstFilters])
jwstFlux[jwstFlux < 1e-90] = 0.0
jwstFlux = np.divide(jwstFlux,jwstNormFlux).tolist()
hubbNormFlux = np.array([np.trapz(hubbleFilters[aFilt](rsWavelnOneAge[::-1])/rsFreqOneAge,rsFreqOneAge).value for aFilt in hubbleFilters])
hubbNormFlux[hubbNormFlux <= 0.0] = -1.0
hubbFlux = np.array([(np.trapz(hubbleFilters[aFilt](rsWavelnOneAge[::-1])*lyForFluxHzOneAge/rsFreqOneAge,rsFreqOneAge)).value
for aFilt in hubbleFilters])
hubbFlux[hubbFlux < 1e-90] = 0.0
hubbFlux = np.divide(hubbFlux,hubbNormFlux).tolist()
jhkNormFlux = np.array([np.trapz(jhkFilters[aFilt](rsWavelnOneAge[::-1])/rsFreqOneAge,rsFreqOneAge).value for aFilt in jhkFilters])
jhkNormFlux[jhkNormFlux <= 0.0] = -1.0
jhkFlux = np.array([(np.trapz(jhkFilters[aFilt](rsWavelnOneAge[::-1])*lyForFluxHzOneAge/rsFreqOneAge,rsFreqOneAge)).value
for aFilt in jhkFilters])
jhkFlux[jhkFlux < 1e-90] = 0.0
jhkFlux = np.divide(jhkFlux,jhkNormFlux).tolist()
aLine = [age] + [z] + jwstFlux + hubbFlux + jhkFlux
outfile.append(aLine)
return outfile
#
# Read in Schaerer and SB99 luminosity files and generate flux in filter.
# These tables only consider Ly-forest absorption -- and redshift.
# User needs to apply reddening.
#
def buildFilterFluxFiles():
"""Creates flux-in-filter files from the SED.
This is untested but a copy of the python notebook that works.
"""
jwstFilters = lf.loadJWSTFilters(suppress=True)
hubbleFilters = lf.loadHubbleFilters(suppress=True)
lamRange = np.logspace(1.95,5.7,5500)
schaererPath = '/Users/earnric/Research/Research-Observability/Software-Models/Schaerer/'
schaererDirs = ['pop3_TA/','pop3_TE/','e-70_mar08/','e-50_mar08/']
Zs = [0.0, 0.0, 1.0e-7, 1.0e-5]
schaererPopFilePattern = 'pop3_ge0_log?_500_001_is5.[0-9]*' # is5 files have ages in step of 1 Myr
schaererLowZFilePattern = 'e-?0_sal_100_001_is2.[0-9]*' # is2 files have ages in step of 0.05 dex
# Load the schaerer files...
# Note that due to spacing (in age), there are sometimes two files that
# are SEDS for the same time-stamp! Skip the second one (and third!)
lastAge = 0.0
for i, (Z, schaererDir) in enumerate(zip(Zs,schaererDirs)):
if schaererDir.startswith('pop3'):
schaererFilePattern = schaererPath + schaererDir + schaererPopFilePattern # Pop III files, 1 Myr spacing
else:
schaererFilePattern = schaererPath + schaererDir + schaererLowZFilePattern # Low Z files, 0.05 dex spacing
schaererFiles = glob.glob(schaererFilePattern) # All the files in the dir...
schaererFiles = [a for a in schaererFiles if not re.search('\.[1-2][0-9][b]*$', a)] # remove .1? and .2? files
schaererAges = [linecache.getline(file,13) for file in schaererFiles] # Get the line with the (log) age...
schaererAges = np.array([float(sa[30:]) for sa in schaererAges],dtype=float) # Log age starts at position 30
schaererData = np.array([np.loadtxt(file,skiprows=16) for file in schaererFiles])
ageSortIndxes = schaererAges.argsort() # Array of indices to sort things by age...
schaererData = schaererData[ageSortIndxes]
schaererAges = schaererAges[ageSortIndxes]
print(len(schaererData))
print(len(schaererAges),schaererAges)
# Ignore data files with the same age! This occurs in the popIII dirs
# because the timestep is smaller than the age-resolution printed in the file
# Hence we get 2 files with different data but the same time stamp
lastAge = 0.0
schaererDataGood = []
schaererAgesGood = []
for ii,(sd,age) in enumerate(zip(schaererData,schaererAges)):
if age == lastAge:
# Remove it
continue
lastAge = age
schaererDataGood.append(sd)
schaererAgesGood.append(float(age))
# The following builds an array of arrays (one for each age) with each array's entries:
# log age, Z, waveln, lum/A
allSchaererData = [np.insert(sed[:,[0,2]],[0],[[anAge] for ii in range(0,len(sed))], axis=1)
for anAge,sed in zip(schaererAgesGood, schaererDataGood)]
allSchaererData = np.array(allSchaererData).reshape(len(allSchaererData)*len(allSchaererData[0]),3)
if i == 0:
pop3TA = allSchaererData # may need a np.copy(...) here... ??
elif i == 1:
pop3TE = allSchaererData
elif i == 2:
Zem7 = allSchaererData
elif i == 3:
Zem5 = allSchaererData
# We now have:
# [[log age, waveln, flux], [], ...]
# Generate the filter-flux data...
redshifts = [2.0,3.0,4.0,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0]
schaererList = [pop3TA, pop3TE, Zem7, Zem5]
schaererNames = ["pop3TA", "pop3TE", "Zem7", "Zem5"]
hdr = 'LogAge, redshift, '
hdr += ', '.join([jFilt for jFilt in jwstFilters])
hdr += ', '
hdr += ', '.join([hFilt for hFilt in hubbleFilters])
hdr += ', '
hdr += ', '.join([jFilt for jFilt in jhkFilters])
fmtStr = '%.4f, %.1f, '
fmtStr += ', '.join(['%.4e' for ii in arange(len(jwstFilters)+len(hubbleFilters)+len(jhkFilters))])
for schaererData, fname in zip(schaererList,schaererNames):
ages = np.unique(schaererData[:,0]) # Get the list of ages
ages = ages[ages <= 9.01] # We don't need flux for stars older than 1.02 Gyr
# Get the log age, wavelength & Luminosity/s/Ang
allAges = schaererData[:,0]
wavelns = schaererData[:,1] * u.Angstrom
LperA = schaererData[:,2] * u.erg / u.second / u.angstrom
for z in redshifts:
buildOutFile(z,allAges,ages,wavelns,LperA)
filename = fname + "_" + str(z) + ".gz"
print('writing {}'.format(filename))
np.savetxt(filename, outfile, fmt=fmtStr, delimiter=', ', header=hdr)
print("Done with Schaerer...")
# Process SB99 files...
# SB99 format: TIME [YR] WAVELENGTH [A] LOG TOTAL LOG STELLAR LOG NEBULAR [ERG/SEC/A]
# REMEMBER, SB99 data is for a population of 1e6 M_sun
SB99Path = '/Users/earnric/OneDrive/STARBURST99-runs/' # Home computer dir...
SB99Dirs = ['padova0004-op/','padova004-op/','padova008-op/','padova02-op/']
Zs = [0.0004, 0.004, 0.008, 0.02]
SB99FilePat = 'padova*.spectrum1'
# SB99 format: TIME [YR] WAVELENGTH [A] LOG TOTAL LOG STELLAR LOG NEBULAR [ERG/SEC/A]
for i, (Z, SB99Dir) in enumerate(zip(Zs,SB99Dirs)):
SB99FilePattern = SB99Path + SB99Dir + SB99FilePat
SB99Files = glob.glob(SB99FilePattern) # All the files in the dir... should be one!
if len(SB99Files) != 1:
print('Error: too many files in an SB99 dir! - ',SB99Path + SB99Dir)
sys.exit()
SB99Data = np.loadtxt(SB99Files[0],skiprows=6)
if i == 0:
SB990004 = np.dstack((np.log10(SB99Data[:,0]),SB99Data[:,1],10**(SB99Data[:,2]-6.0))).reshape(len(SB99Data[:,1]),3)
elif i == 1:
SB99004 = np.dstack((np.log10(SB99Data[:,0]),SB99Data[:,1],10**(SB99Data[:,2]-6.0))).reshape(len(SB99Data[:,1]),3)
elif i == 2:
SB99008 = np.dstack((np.log10(SB99Data[:,0]),SB99Data[:,1],10**(SB99Data[:,2]-6.0))).reshape(len(SB99Data[:,1]),3)
elif i == 3:
SB9902 = np.dstack((np.log10(SB99Data[:,0]),SB99Data[:,1],10**(SB99Data[:,2]-6.0))).reshape(len(SB99Data[:,1]),3)
# We now have:
# [[log age, waveln, flux], [], ...]
# Generate the filter-flux data...
redshifts = [2.0,3.0,4.0,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0]
# redshifts = [14.0,15.0,16.0]
sb99List = [SB990004, SB99004, SB99008, SB9902]
sb99Names = ["SB990004", "SB99004", "SB99008", "SB9902"]
# sb99List = [SB9902]
# sb99Names = ["SB9902"]
print(hdr)
for sb99Data, fname in zip(sb99List,sb99Names):
ages = np.unique(sb99Data[:,0]) # Get the list of ages
# Get the log age, wavelength & Luminosity/s/Ang
allAges = sb99Data[:,0]
wavelns = sb99Data[:,1] * u.Angstrom
LperA = sb99Data[:,2] * u.erg / u.second / u.angstrom
for z in redshifts:
buildOutFile(z,allAges,ages,wavelns,LperA)
filename = fname + "_" + str(z) + ".gz"
print('writing {}'.format(filename))
np.savetxt(filename, outfile, fmt=fmtStr, delimiter=', ', header=hdr)
print("Done with SB99...")
return