-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmake_sim_suite.py
377 lines (267 loc) · 11.7 KB
/
make_sim_suite.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
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
import time as t
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from pixell import enmap, reproject, utils, curvedsky, wcsutils, bunch
from orphics import maps, io
import healpy as hp
from past.utils import old_div
import argparse
import sys
start = t.time()
paths = bunch.Bunch(io.config_from_yaml("input/sim_data.yml"))
print("Paths: ", paths)
# RUN SETTING ------------------------------------------------------------------
parser = argparse.ArgumentParser()
parser.add_argument(
"save_name", type=str, help="Name you want for your output."
)
parser.add_argument(
"which_sim", type=str, help="Choose the sim e.g. websky or sehgal or agora."
)
parser.add_argument(
"--high-accuracy", action="store_true", help="If using a high accuracy websky map, choose this option."
)
args = parser.parse_args()
output_path = paths.simsuite_path
save_name = args.save_name
if args.which_sim == "websky":
print(" ::: producing maps for WEBSKY sim")
save_name = "websky_" + save_name
sim_path = paths.websky_sim_path
elif args.which_sim == "sehgal":
print(" ::: producing maps for SEHGAL sim")
save_name = "sehgal_" + save_name
sim_path = paths.sehgal_sim_path
elif args.which_sim == "agora":
print(" ::: producing maps for AGORA sim")
save_name = "agora_" + save_name
sim_path = paths.agora_sim_path
save_dir = f"{output_path}/{save_name}/"
io.mkdir(f"{save_dir}")
print(" ::: saving to", save_dir)
# SIM SETTING ------------------------------------------------------------------
freq_sz = 150
print(" ::: map frequency at %d GHz" %freq_sz)
px = 0.5
fwhm = 1.5
nlevel = 15.0
# BEAM CONVOLUTION -------------------------------------------------------------
def apply_beam(imap):
# map2alm of the maps, almxfl(alm, beam_1d) to convolve with beam, alm2map to convert back to map
if args.which_sim == "websky" or args.which_sim == "agora": nside = 8192
elif args.which_sim == "sehgal": nside = 4096
alm_lmax = nside * 3
bfunc = lambda x: maps.gauss_beam(fwhm, x)
imap_alm = curvedsky.map2alm(imap, lmax=alm_lmax)
beam_convoloved_alm = curvedsky.almxfl(imap_alm, bfunc)
return curvedsky.alm2map(beam_convoloved_alm, enmap.empty(imap.shape, imap.wcs))
# PREPARING MAPS ---------------------------------------------------------------
# paths for reprojected maps
if args.which_sim == "websky":
if args.high_accuracy:
r_lmap = sim_path + paths.websky_dlensed_reproj
r_kmap = sim_path + paths.websky_kappa4p5_reproj
r_tszmap = sim_path + paths.websky_dtsz_reproj
r_kszmap = sim_path + paths.websky_dksz_reproj
r_cib093map = sim_path + paths.websky_dcib093_reproj
r_cib145map = sim_path + paths.websky_dcib145_reproj
else:
r_lmap = sim_path + paths.websky_cmb_reproj
r_kmap = sim_path + paths.websky_kappa_reproj
r_tszmap = sim_path + paths.websky_tsz_reproj
elif args.which_sim == "sehgal":
r_lmap = sim_path + paths.sehgal_cmb_reproj
r_kmap = sim_path + paths.sehgal_kappa_reproj
r_tszmap = sim_path + paths.sehgal_tsz_reproj
elif args.which_sim == "agora":
r_lmap = sim_path + paths.agora_cmb_reproj
r_kmap = sim_path + paths.agora_kappa_reproj
r_tszmap = sim_path + paths.agora_tsz_reproj
# reading lensed cmb map
try:
lmap = enmap.read_map(r_lmap, delayed=False)
print(" ::: reading reprojected lensed cmb map:", r_lmap)
except:
shape, wcs = enmap.fullsky_geometry(res=px*utils.arcmin, proj="car")
if args.which_sim == "websky":
if args.high_accuracy:
ifile = f"{paths.mat_path}/dlensed.fits"
lmap = enmap.read_map(ifile)
print(" ::: reading lensed cmb map:", ifile)
else:
ifile = f"{sim_path}lensed_alm.fits"
alm = np.complex128(hp.read_alm(ifile, hdu=(1, 2, 3)))
lmap = curvedsky.alm2map(alm[0,:], enmap.empty(shape, wcs, dtype=np.float64))
print(" ::: reading lensed alm map and converting to lensed cmb map:", ifile)
elif args.which_sim == "sehgal":
ifile = f"{sim_path}Sehgalsimparams_healpix_4096_KappaeffLSStoCMBfullsky_phi_SimLens_Tsynfastnopell_fast_lmax8000_nside4096_interp2.5_method1_1_lensed_map.fits"
hmap = hp.read_map(ifile).astype(np.float64)
lmap = reproject.healpix2map(hmap, shape, wcs)
print(" ::: reading lensed cmb map:", ifile)
elif args.which_sim == "agora":
ifile = f"{sim_path}cmb/len/tqu1/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap.fits"
hmap,_,_ = hp.read_map(ifile, field=[0,1,2]).astype(np.float64)
lmap = reproject.healpix2map(hmap, shape, wcs)
print(" ::: reading lensed cmb map:", ifile)
enmap.write_map(r_lmap, lmap)
print("lmap", lmap.shape)
shape, wcs = lmap.shape, lmap.wcs
# reading true kappa map
try:
kmap = enmap.read_map(r_kmap, delayed=False)
print(" ::: reading reprojected true kappa map:", r_kmap)
except:
if args.which_sim == "websky":
if args.high_accuracy:
ifile = f"{sim_path}kap_lt4.5.fits" # CMB lensing convergence from z<4.5 from halo+field websky
else:
ifile = f"{sim_path}kap.fits" # CMB lensing convergence from 0<z<1100
hmap = np.atleast_2d(hp.read_map(ifile, field=tuple(range(0,1)))).astype(np.float64)
kmap = reproject.healpix2map(hmap, shape, wcs)[0,:,:]
elif args.which_sim == "sehgal":
ifile = f"{sim_path}healpix_4096_KappaeffLSStoCMBfullsky.fits"
hmap = hp.read_map(ifile).astype(np.float64)
kmap = reproject.healpix2map(hmap, shape, wcs)
elif args.which_sim == "agora":
ifile = f'{sim_path}../../kappa_alm_recon_agora_phiNG_phi1_seed1.fits'
hmap = hp.read_alm(ifile).astype(np.complex128)
kmap = curvedsky.alm2map(hmap, enmap.empty(shape, wcs, dtype=np.float64))
enmap.write_map(r_kmap, kmap)
print(" ::: reading and reprojecting true kappa map:", ifile)
print("kmap", kmap.shape)
# reading tsz map
try:
tszmap = enmap.read_map(r_tszmap, delayed=False)
print(" ::: reading reprojected tsz map:", r_tszmap)
except:
if args.which_sim == "websky":
ifile = f"{sim_path}tsz_8192.fits"
hmap = np.atleast_2d(hp.read_map(ifile, field=tuple(range(0,1)))).astype(np.float64)
ymap = reproject.healpix2map(hmap, shape, wcs)[0,:,:]
elif args.which_sim == "sehgal":
ifile = f"{sim_path}tSZ_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75.fits"
hmap = hp.read_map(ifile).astype(np.float64)
ymap = reproject.healpix2map(hmap, shape, wcs)
elif args.which_sim == "agora":
ifile = f"{sim_path}tsz/len/agora_ltszNG_bahamas78_bnd_unb_1.0e+12_1.0e+18_lensed.fits"
hmap = hp.read_map(ifile).astype(np.float64)
ymap = reproject.healpix2map(hmap, shape, wcs)
print(" ::: reading ymap:", ifile)
print("ymap", ymap.shape)
# convert compton-y to delta-T (in uK)
tcmb = 2.726
tcmb_uK = tcmb * 1e6 #micro-Kelvin
H_cgs = 6.62608e-27
K_cgs = 1.3806488e-16
def fnu(nu):
"""
nu in GHz
tcmb in Kelvin
"""
mu = H_cgs*(1e9*nu)/(K_cgs*tcmb)
ans = mu/np.tanh(old_div(mu,2.0)) - 4.0
return ans
tszmap = fnu(freq_sz) * ymap * tcmb_uK
print(" ::: converting ymap to tsz map at %d GHz" %freq_sz)
enmap.write_map(r_tszmap, tszmap)
print("tszmap", tszmap.shape)
# reading ksz map
try:
kszmap = enmap.read_map(r_kszmap, delayed=False)
print(" ::: reading reprojected ksz map:", r_kszmap)
except:
if args.which_sim == "websky":
ifile = f"{sim_path}ksz.fits"
hmap = np.atleast_2d(hp.read_map(ifile, field=tuple(range(0,1)))).astype(np.float64)
kszmap = reproject.healpix2map(hmap, shape, wcs)[0,:,:]
enmap.write_map(r_kszmap, kszmap)
print(" ::: reading and reprojecting ksz map:", ifile)
print("kszmap", kszmap.shape)
# reading cib map
try:
if freq_sz == 90: r_cibmap = r_cib093map
elif freq_sz == 150: r_cibmap = r_cib145map
cibmap = enmap.read_map(r_cibmap, delayed=False)
print(" ::: reading reprojected cib map:", r_cibmap)
except:
if args.which_sim == "websky":
# frequency for CIB
if freq_sz == 90: ifile = f"{sim_path}cib_nu0093.fits"
elif freq_sz == 150: ifile = f"{sim_path}cib_nu0145.fits"
hmap = np.atleast_2d(hp.read_map(ifile, field=tuple(range(0,1)))).astype(np.float64)
cibmap_i = reproject.healpix2map(hmap, shape, wcs)[0,:,:]
print(" ::: reading cib map:", ifile)
print("cibmap", np.shape(cibmap_i))
kboltz = 1.3806503e-23 #MKS
hplanck = 6.626068e-34 #MKS
clight = 299792458.0 #MKS
def ItoDeltaT(nu):
# conversion from specific intensity to Delta T units (i.e., 1/dBdT|T_CMB)
# i.e., from W/m^2/Hz/sr (1e-26 Jy/sr) --> uK_CMB
# i.e., you would multiply a map in 1e-26 Jy/sr by this factor to get an output map in uK_CMB
nu *= 1e9
X = hplanck*nu/(kboltz*tcmb)
dBnudT = (2.*hplanck*nu**3.)/clight**2. * (np.exp(X))/(np.exp(X)-1.)**2. * X/tcmb_uK * 1e26
return 1./dBnudT
# frequency for CIB
if freq_sz == 150: freq_cib = 145
elif freq_sz == 90: freq_cib = 93
cibmap = ItoDeltaT(freq_cib) * cibmap_i
print(" ::: converting cib map at %d GHz" %freq_cib)
enmap.write_map(r_cibmap, cibmap)
print("cibmap", cibmap.shape)
assert wcsutils.equal(kmap.wcs, lmap.wcs)
assert wcsutils.equal(kmap.wcs, tszmap.wcs)
assert wcsutils.equal(kmap.wcs, kszmap.wcs)
assert wcsutils.equal(kmap.wcs, cibmap.wcs)
scmb = lmap
scmb_cib = lmap + cibmap
scmb_ksz_cib = lmap + kszmap + cibmap
scmb_tsz = lmap + tszmap
scmb_tsz_cib = lmap + tszmap + cibmap
scmb_tsz_ksz_cib = lmap + tszmap + kszmap + cibmap
print(" ::: SIGNAL maps are ready!")
white_noise = maps.white_noise(lmap.shape, lmap.wcs, noise_muK_arcmin=nlevel)
print(" ::: applying beam and adding white noise of %.1f uK" %nlevel)
bcmb = apply_beam(scmb)
bcmb_ksz_cib = apply_beam(scmb_ksz_cib)
ocmb = apply_beam(scmb) + white_noise
ocmb_cib = apply_beam(scmb_cib) + white_noise
ocmb_ksz_cib = apply_beam(scmb_ksz_cib) + white_noise
ocmb_tsz = apply_beam(scmb_tsz) + white_noise
ocmb_tsz_cib = apply_beam(scmb_tsz_cib) + white_noise
ocmb_tsz_ksz_cib = apply_beam(scmb_tsz_ksz_cib) + white_noise
print(" ::: OBSERVED maps are ready!")
save_scmb = f"{save_dir}scmb.fits"
save_scmb_cib = f"{save_dir}scmb_cib.fits"
save_scmb_ksz_cib = f"{save_dir}scmb_ksz_cib.fits"
save_scmb_tsz = f"{save_dir}scmb_tsz.fits"
save_scmb_tsz_cib = f"{save_dir}scmb_tsz_cib.fits"
save_scmb_tsz_ksz_cib = f"{save_dir}scmb_tsz_ksz_cib.fits"
save_ocmb = f"{save_dir}ocmb.fits"
save_ocmb_cib = f"{save_dir}ocmb_cib.fits"
save_ocmb_ksz_cib = f"{save_dir}ocmb_ksz_cib.fits"
save_ocmb_tsz = f"{save_dir}ocmb_tsz.fits"
save_ocmb_tsz_cib = f"{save_dir}ocmb_tsz_cib.fits"
save_ocmb_tsz_ksz_cib = f"{save_dir}ocmb_tsz_ksz_cib.fits"
save_bcmb = f"{save_dir}bcmb.fits"
save_bcmb_ksz_cib = f"{save_dir}bcmb_ksz_cib.fits"
print(" ::: saving all maps")
enmap.write_map(save_scmb, scmb)
enmap.write_map(save_scmb_cib, scmb_cib)
enmap.write_map(save_scmb_ksz_cib, scmb_ksz_cib)
enmap.write_map(save_scmb_tsz, scmb_tsz)
enmap.write_map(save_scmb_tsz_cib, scmb_tsz_cib)
enmap.write_map(save_scmb_tsz_ksz_cib, scmb_tsz_ksz_cib)
enmap.write_map(save_ocmb, ocmb)
enmap.write_map(save_ocmb_cib, ocmb_cib)
enmap.write_map(save_ocmb_ksz_cib, ocmb_ksz_cib)
enmap.write_map(save_ocmb_tsz, ocmb_tsz)
enmap.write_map(save_ocmb_tsz_cib, ocmb_tsz_cib)
enmap.write_map(save_ocmb_tsz_ksz_cib, ocmb_tsz_ksz_cib)
enmap.write_map(save_bcmb, bcmb)
enmap.write_map(save_bcmb_ksz_cib, bcmb_ksz_cib)
elapsed = t.time() - start
print("\r ::: entire run took %.1f seconds" %elapsed)