-
-
Save erussier/7527cde4caf32887dd24379573b4d923 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import healpy as hp | |
import matplotlib.pyplot as plt | |
import numpy as np | |
import pymaster as nmt | |
from astropy.io import fits | |
import pysm3 as ps | |
import pysm3.units as u | |
nside = 2048 | |
hpx_datapath = './Healpix_3.82/data/' | |
npix = hp.nside2npix(nside) | |
band = 217 | |
res_217 = 4.66 | |
#Smooth map to 30 arcmin | |
def resol_change(data_path, nside, old_res_ang_arcmin, new_res_ang_deg, hpx_datapath, conv_data_map): | |
lmax = 3*nside - 1 | |
data_map = hp.read_map(data_path, field=(0, 1, 2)) | |
alms = hp.map2alm(data_map, lmax=lmax, use_weights=True, datapath=hpx_datapath) | |
Bl_data_map = hp.gauss_beam(np.deg2rad(old_res_ang_arcmin/60.), lmax=lmax, pol=True) | |
Bl_com = hp.gauss_beam(np.deg2rad(new_res_ang_deg), lmax=lmax, pol=True) | |
alms[0] = hp.almxfl(alms[0], Bl_com[:,0]/Bl_data_map[:,0]) | |
alms[1] = hp.almxfl(alms[1], Bl_com[:,1]/Bl_data_map[:,1]) | |
alms[2] = hp.almxfl(alms[2], Bl_com[:,2]/Bl_data_map[:,2]) | |
data_map = hp.alm2map(alms, nside) | |
data_map = data_map * conv_data_map | |
return data_map | |
#Models | |
d9_bpi_path = f"./data/dust9_HFI{band}-bpi_wo_renorm_" + str(nside) + '.fits' | |
d10_bpi_path = f"./data/dust10_HFI{band}-bpi_wo_renorm_" + str(nside) + '.fits' | |
d9_bpi_map = hp.read_map(d9_bpi_path, field = None) | |
d10_bpi_map = hp.read_map(d10_bpi_path, field = None) | |
I_d9_bpi = d9_bpi_map[0] | |
I_d10_bpi = d10_bpi_map[0] | |
Ipol_d9_bpi = np.sqrt(d9_bpi_map[1]**2 + d9_bpi_map[2]**2) | |
Ipol_d10_bpi = np.sqrt(d10_bpi_map[1]**2 + d10_bpi_map[2]**2) | |
d9_bpi_map_res_change = resol_change(d9_bpi_path, nside, res_217, 0.5, hpx_datapath, 1) | |
print("end d9 res change") | |
d10_bpi_map_res_change = resol_change(d10_bpi_path, nside, res_217, 0.5, hpx_datapath, 1) | |
print("end d10 res change") | |
#Save maps | |
c1 = fits.Column(name='I', array=d9_bpi_map_res_change[0], format='D') | |
c2 = fits.Column(name='Q', array=d9_bpi_map_res_change[1], format='D') | |
c3 = fits.Column(name='U', array=d9_bpi_map_res_change[2], format='D') | |
t = fits.BinTableHDU.from_columns([c1, c2, c3]) | |
t.writeto(f'./data/{band}_d9_bpi_res_30arcmin_2048.fits', overwrite = True) | |
print("end write file d9 res change") | |
c1 = fits.Column(name='I', array=d10_bpi_map_res_change[0], format='D') | |
c2 = fits.Column(name='Q', array=d10_bpi_map_res_change[1], format='D') | |
c3 = fits.Column(name='U', array=d10_bpi_map_res_change[2], format='D') | |
t = fits.BinTableHDU.from_columns([c1, c2, c3]) | |
t.writeto(f'./data/{band}_d10_bpi_res_30arcmin_2048.fits', overwrite = True) | |
print("end write file d10 res change") | |
#Data | |
# pr3_hfi_path = "/global/cfs/cdirs/cmb/data/planck2018/pr3/frequencymaps/HFI_SkyMap_217_2048_R3.01_full.fits" | |
# noise_sim_path = "/global/cfs/cdirs/cmb/data/planck2018/ffp10/mc_noise/217/ffp10_noise_217_full_map_mc_00000.fits" | |
# pr3_hfi_353_map = hp.read_map(pr3_hfi_path, field = None)*1e6 | |
# I_pr3_hfi_353 = pr3_hfi_353_map[0] | |
# Ipol_pr3_hfi_353 = np.sqrt(pr3_hfi_353_map[1]**2 + pr3_hfi_353_map[2]**2) | |
# pr3_res_change = resol_change(pr3_hfi_path, nside, res_217, 0.5, hpx_datapath, 1e6) | |
# print("end pr3 res change") | |
# c1 = fits.Column(name='I', array=pr3_res_change[0], format='D') | |
# c2 = fits.Column(name='Q', array=pr3_res_change[1], format='D') | |
# c3 = fits.Column(name='U', array=pr3_res_change[2], format='D') | |
# t = fits.BinTableHDU.from_columns([c1, c2, c3]) | |
# t.writeto(f'./data/{band}_pr3_res_30arcmin_2048.fits', overwrite = True) | |
# print("end write file pr3 res change") | |
# noise_sim_353, hd_noise_sim_353 = hp.read_map(noise_sim_path, h = True, field = None) | |
# noise_sim_uKcmb = noise_sim_353*1e6 | |
# noise_res_change = resol_change(noise_sim_path, nside, res_857, 0.5, hpx_datapath, 1e6) | |
# print("end noise res change") | |
# c1 = fits.Column(name='I', array=noise_res_change[0], format='D') | |
# c2 = fits.Column(name='Q', array=noise_res_change[1], format='D') | |
# c3 = fits.Column(name='U', array=noise_res_change[2], format='D') | |
# t = fits.BinTableHDU.from_columns([c1, c2, c3]) | |
# t.writeto(f'./data/{band}_noise_res_30arcmin_2048.fits', overwrite = True) | |
# print("end write file noise res change") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment