857_comp_dust_map_models.py

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 = 857
res_857 = 4.23

#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=None)
    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 = "./data/dust9_HFI857-bpi_wo_renorm_" + str(nside) + '.fits'
d10_bpi_path = "./data/dust10_HFI857-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)
#30 arcmin maps
d9_bpi_map_res_change = resol_change(d9_bpi_path, nside, res_857, 0.5, hpx_datapath, 1)
print("end d9 res change")
d10_bpi_map_res_change = resol_change(d10_bpi_path, nside, res_857, 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_857_2048_R3.01_full.fits"
# # noise_sim_path = "/global/cfs/cdirs/cmb/data/planck2018/ffp10/mc_noise/857/ffp10_noise_857_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_857, 0.5, hpx_datapath, 1)
# 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")