bandpass_int_dust.py

import warnings
warnings.filterwarnings("ignore") 

import pysm3 as pysm
import pysm3.units as u
import numpy as np
import healpy as hp 
import matplotlib.pyplot as plt
import scipy.stats as stat
import matplotlib as mpl
import astropy.io.fits as fits
import scipy.constants as con

mpl.rcParams['font.sans-serif'] = "Helvetica"
plt.style.use('default')

nside = 2048
npix = hp.nside2npix(nside)
lmax = 3 * nside - 1 
hpx_datapath = './Healpix_3.82/data/'

RIMO_file = './data/HFI_RIMO_R3.00.fits'
hdul = fits.open(RIMO_file)

band = '353'

nus_in_GHz = hdul['BANDPASS_F'+band].data['WAVENUMBER']*con.c*1e-7
transmission = hdul['BANDPASS_F'+band].data['TRANSMISSION']

nus_in_GHz = nus_in_GHz[transmission>0.001*np.max(transmission)]
transmission = transmission[transmission>0.001*np.max(transmission)]
transmission /= np.trapz(transmission, x=nus_in_GHz*con.giga)
print("end_read_file + renormalization", "\n")
# transmission /= np.trapz(transmission, x=nus_in_GHz*con.giga)
nus_in_GHz = np.array(nus_in_GHz, dtype=np.float64)
transmission = np.array(transmission,dtype=np.float64)
print("end create array", "\n")
UC = pysm.bandpass_unit_conversion(nus_in_GHz*u.GHz, weights=transmission, output_unit=u.uK_CMB, input_unit=u.uK_RJ, cut=1e-10)
print("UC finished", "\n")

dust_12 = pysm.Sky(nside=nside, preset_strings=["d12"])
dustmap_12 = dust_12.get_emission(nus_in_GHz * u.GHz, transmission)
print("d12 get emission finished", "\n")
dustmap_12 = dustmap_12 * UC
dustmap_12 = hp.smoothing(dustmap_12, fwhm=np.deg2rad(4.41/60.), use_weights=True, datapath=hpx_datapath)
hp.write_map('./data/dust12_HFI353-bpi_renorm_' + str(nside) + '.fits', dustmap_12, dtype=['float32','float32','float32'], overwrite=True)
print("d12 write file finished", "\n")
dust_9 = pysm.Sky(nside=nside, preset_strings=["d9"])
dustmap_9 = dust_9.get_emission(nus_in_GHz * u.GHz, transmission)
print("d9 get emission finished", "\n")
dustmap_9 = dustmap_9 * UC
dustmap_9 = hp.smoothing(dustmap_9, fwhm=np.deg2rad(4.41/60.), use_weights=True, datapath=hpx_datapath)
hp.write_map('./data/dust9_HFI353-bpi_renorm_' + str(nside) + '.fits', dustmap_9, dtype=['float32','float32','float32'], overwrite=True)
print("d9 write file finished", "\n")
dust_10 = pysm.Sky(nside=nside, preset_strings=["d10"])
dustmap_10 = dust_10.get_emission(nus_in_GHz * u.GHz, transmission)
print("d10 get emission finished", "\n")
dustmap_10 = dustmap_10 * UC
dustmap_10 = hp.smoothing(dustmap_10, fwhm=np.deg2rad(4.41/60.), use_weights=True, datapath=hpx_datapath)
hp.write_map('./data/dust10_HFI353-bpi_renorm_' + str(nside) + '.fits', dustmap_10, dtype=['float32','float32','float32'], overwrite=True)
print("d10 write file finished", "\n")
dust_11 = pysm.Sky(nside=nside, preset_strings=["d11"])
dustmap_11 = dust_11.get_emission(nus_in_GHz * u.GHz, transmission)
print("d11 get emission finished", "\n")
dustmap_11 = dustmap_11 * UC
dustmap_11 = hp.smoothing(dustmap_11, fwhm=np.deg2rad(4.41/60.), use_weights=True, datapath=hpx_datapath)
hp.write_map('./data/dust11_HFI353-bpi_renorm_' + str(nside) + '.fits', dustmap_11, dtype=['float32','float32','float32'], overwrite=True)
print("d11 write file finished", "\n")