-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathinput_var.py
232 lines (199 loc) · 11 KB
/
input_var.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
from headers_constants import *
class data_var(object):
def __init__(self, exp, mass, z, ell):
# ############### cib data #########################
self.exp = exp
name = self.exp['name']
deltah_cib = 200
self.z_c = 1.5
self.freqcib = self.exp['freq_cib']
self.cc = self.exp['cc']
self.fc = self.exp['fc']
self.cc_cibmean = self.exp['cc_cibmean']
self.freq_cibmean = self.exp['freq_cibmean']
self.mass = mass
self.z = z
self.ell = ell
nm = len(self.mass)
nz = len(self.z)
# ########## reading in the matter power spectra #############
redshifts = np.loadtxt('data_files/redshifts.txt')
if min(self.z) < min(redshifts) or max(self.z) > max(redshifts):
print ("If the redshift range is outside of [%s to %s], then " +
"values of the matter power spectrum and effective " +
"CIB SEDs are extrapolated and might be incorrect.") % (min(redshifts), max(redshifts))
ll = [str(x) for x in range(1, 211)]
"""
please note that the matter power spectrum files here are arranged in
a reversed order i.e. redshift decreases as you go from _1 file to _210.
Perhaps better to generate your own power spectra for the redshift
array you want to consider and read that here. Also, note that I am
getting rid of the reduced Hubble constant here from units of k and Pk.
"""
addr = 'data_files/matter_power_spectra'
pkarray = np.loadtxt('%s/test_highk_lin_matterpower_210.dat' % (addr))
k = pkarray[:, 0]*cosmo.h
Pk = np.zeros((len(k), len(redshifts)))
for i in range(len(redshifts)):
pkarray = np.loadtxt("%s/test_highk_lin_matterpower_%s.dat" % (addr, ll[209-i]))
Pk[:, i] = pkarray[:, 1]/cosmo.h**3
pkinterpz = interp1d(redshifts, Pk, kind='linear', bounds_error=False, fill_value="extrapolate")
self.k_array = np.zeros((len(self.ell), len(self.z)))
self.Pk_int = np.zeros(self.k_array.shape)
"""
Pk_int 2-d array for corresponding redshifts and
given ell range such that k = ell/chi i.e. for every
redshift
"""
chiz = cosmo.comoving_distance(self.z).value
for i in range(len(self.ell)):
self.k_array[i, :] = self.ell[i]/chiz
for j in range(len(self.z)):
pkz = pkinterpz(self.z[j])
self.Pk_int[i, j] = np.interp(self.k_array[i, j], k, pkz)
if self.exp['do_cib'] == 1 or self.exp['do_cibxtsz'] == 1:
# ######### reading and interpolating the SEDs
"""
The effective SEDs for the CIB for Planck (100, 143, 217, 353, 545,
857) and
IRAS (3000) GHz frequencies. Taking out the IRAS SEDs and that is
why snu[:-1, :] in the last line here.
Here we are shwoing the CIB power spectra corressponding to the
Planck
frequency channels. If you want to calculate the Hershel/Spire
power spectra, use corresponding files in the data folder.
"""
if name == 'Planck':
snuaddr = 'data_files/filtered_snu_planck.fits'
hdulist = fits.open(snuaddr)
redshifts = hdulist[1].data
snu_eff = hdulist[0].data # in Jy/Lsun
hdulist.close()
fsnu_eff = interp1d(redshifts, snu_eff, kind='linear',
bounds_error=False, fill_value="extrapolate")
self.snu = fsnu_eff(self.z)[:-1, :]
elif name == 'Herschel-spire':
snuaddr = 'data_files/filtered_snu_spire.fits'
hdulist = fits.open(snuaddr)
redshifts = hdulist[1].data
snu_eff = hdulist[0].data # in Jy/Lsun
hdulist.close()
fsnu_eff = interp1d(redshifts, snu_eff, kind='linear',
bounds_error=False, fill_value="extrapolate")
self.snu = fsnu_eff(self.z)
else:
# ######### unfiltered SEDs ###########################
list_of_files = sorted(glob.glob('data_files/TXT_TABLES_2015/./*.txt'))
a = list_of_files[95]
b = list_of_files[96]
for i in range(95, 208):
list_of_files[i] = list_of_files[i+2]
list_of_files[208] = a
list_of_files[209] = b
wavelengths = np.loadtxt('data_files/TXT_TABLES_2015/EffectiveSED_B15_z0.012.txt')[:, [0]]
# wavelengths are in microns
freq = c_light/wavelengths
# c_light is in Km/s, wavelength is in microns and we would like to
# have frequency in GHz. So gotta multiply by the following
# numerical factor which comes out to be 1
# numerical_fac = 1e3*1e6/1e9
numerical_fac = 1.
freqhz = freq*1e3*1e6
freq *= numerical_fac
freq_rest = freqhz*(1+redshifts)
n = np.size(wavelengths)
snu_unfiltered = np.zeros([n, len(redshifts)])
for i in range(len(list_of_files)):
snu_unfiltered[:, i] = np.loadtxt(list_of_files[i])[:, 1]
L_IR15 = self.L_IR(snu_unfiltered, freq_rest, redshifts)
# print (L_IR15)
for i in range(len(list_of_files)):
snu_unfiltered[:, i] = snu_unfiltered[:, i]*L_sun/L_IR15[i]
# Currently unfiltered snus are ordered in increasing wavelengths,
# we re-arrange them in increasing frequencies i.e. invert it
freq = freq[::-1]
snu_unfiltered = snu_unfiltered[::-1]
fsnu_unfiltered = RectBivariateSpline(freq, redshifts,
snu_unfiltered)
nuinp = self.freqcib
self.snu = fsnu_unfiltered(nuinp, self.z)
if self.exp['do_cib'] == 1:
# ######### CIB halo model parameters ###################
cibparresaddr = 'data_files/one_halo_bestfit_allcomponents_lognormal_sigevol_1p5zcutoff_nospire_fcpl_onlyautoshotpar_no3000_gaussian600n857n1200_planck_spire_hmflog10.txt'
self.Meffmax, self.etamax, self.sigmaMh, self.tau = np.loadtxt(cibparresaddr)[:4, 0]
# self.Meffmax, self.etamax, self.sigmaMh, self.tau = 8753289339381.791, 0.4028353504978569, 1.807080723258688, 1.2040244128818796
# ######## hmf, bias, nfw ###########
print ("Calculating the halo mass function, halo bias, nfw " +
"profile " +
"for given mass and redshift for CIB calculations.")
self.hmf_cib = np.zeros((nm, nz))
self.u_nfw_cib = np.zeros((nm, len(self.k_array[:, 0]), nz))
self.bias_m_z_cib = np.zeros((nm, nz))
delta_h = deltah_cib
for r in range(nz):
pkz = pkinterpz(self.z[r])
instance = hmf_unfw_bias.h_u_b(k, pkz, self.z[r],
cosmo, delta_h, self.mass)
self.hmf_cib[:, r] = instance.dn_dlogm()
# nfw_u[:, :, r] = instance.nfwfourier_u()
self.bias_m_z_cib[:, r] = instance.b_nu()
instance2 = hmf_unfw_bias.h_u_b(self.k_array[:, r],
self.Pk_int[:, r], self.z[r],
cosmo, delta_h, self.mass)
self.u_nfw_cib[:, :, r] = instance2.nfwfourier_u()
if self.exp['do_tsz'] == 1 or self.exp['do_cibxtsz'] == 1:
# ############################### tSZ params #####################
xstep = 50
lnx = np.linspace(-6, 1, xstep)
self.x = 10**lnx
# self.nutsz = np.array([100., 143., 217., 353., 545., 857.])*ghz
self.nutsz = np.array(self.freqcib)*ghz
#nus = ['100', '143', '217', '353', '545', '857']
self.delta_h_tsz = 500 # 500 # 200
self.B = 1.5 # 1.41
self.m500 = np.repeat(self.mass[..., np.newaxis], len(self.z),
axis=1)
print ("Calculating the halo mass function, halo bias, nfw " +
"profile " +
"for given mass and redshift for tSZ calculations.")
self.hmf_tsz = np.zeros((len(self.m500), nz))
self.bias_m_z_tsz = np.zeros((len(self.m500[:, 0]), nz))
self.u_nfw_tsz = np.zeros((nm, len(self.k_array[:, 0]), nz))
delta_h = self.delta_h_tsz
for r in range(nz):
pkz = pkinterpz(self.z[r])
instance = hmf_unfw_bias.h_u_b(k, pkz, self.z[r],
cosmo, delta_h, self.m500[:, 0])
self.hmf_tsz[:, r] = instance.dn_dlogm()
# nfw_u[:, :, r] = instance.nfwfourier_u()
self.bias_m_z_tsz[:, r] = instance.b_nu()
instance2 = hmf_unfw_bias.h_u_b(self.k_array[:, r],
self.Pk_int[:, r], self.z[r],
cosmo, delta_h, self.mass)
self.u_nfw_tsz[:, :, r] = instance2.nfwfourier_u()
if self.exp['do_cibxtsz'] == 1:
# mass definition here is m500 for both the CIB and tSZ. The best fit values
# for CIB parameters change due to this.
# self.snu = self.snu #[:-1, :]
# self.cc = self.exp['cc']#[:-1]
# self.fc = self.exp['fc']#[:-1]
# ################################ cib x tSZ #####################
cibxtszparresaddr = 'data_files/one_halo_bestfit_allcomponents_lognormal_sigevol_highk_deltah500_onlyautoshotpar_no3000_gaussian600n857n1200_planck_spire_hmflog10.txt'
# self.Meffmax, self.etamax, self.sigmaMh, self.tau = np.loadtxt(cibxtszparresaddr)[:4, 0]
# self.Meffmax_cross, self.etamax_cross, self.sigmaMh_cross = 6962523672799.227, 0.4967291547804018, 1.8074450009861387
# self.tau_cross = 1.2016980179374213
def L_IR(self, snu_eff, freq_rest, redshifts):
# freq_rest *= ghz # GHz to Hz
fmax = 3.7474057250000e13 # 8 micros in Hz
fmin = 2.99792458000e11 # 1000 microns in Hz
no = 10000
fint = np.linspace(np.log10(fmin), np.log10(fmax), no)
L_IR_eff = np.zeros((len(redshifts)))
dfeq = np.array([0.]*no, dtype=float)
for i in range(len(redshifts)):
L_feq = snu_eff[:, i]*4*np.pi*(Mpc_to_m*cosmo.luminosity_distance(redshifts[i]).value)**2/(w_jy*(1+redshifts[i]))
Lint = np.interp(fint, np.log10(np.sort(freq_rest[:, i])),
L_feq[::-1])
dfeq = 10**(fint)
L_IR_eff[i] = np.trapz(Lint, dfeq)
return L_IR_eff