bandpass_int_dust_857.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 = '857'
res_band = 4.23#arcmin for 857
#res_band = 4.66#arcmin for 217

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")
#for 857GHz 
UC = pysm.bandpass_unit_conversion(nus_in_GHz*u.GHz, weights=transmission, output_unit = u.MJy/ u.sr, input_unit=u.uK_RJ, cut=1e-10)
#for 353, 217GHz
#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(res_band/60.), use_weights=True, datapath=hpx_datapath)
# hp.write_map(f'./data/dust12_HFI{band}-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(res_band/60.), use_weights=True, datapath=hpx_datapath)
# hp.write_map(f'./data/dust9_HFI{band}-bpi_wo_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(res_band/60.), use_weights=True, datapath=hpx_datapath)
hp.write_map(f'./data/dust10_HFI{band}-bpi_wo_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(res_band/60.), use_weights=True, datapath=hpx_datapath)
# hp.write_map(f'./data/dust11_HFI{band}-bpi_renorm_' + str(nside) + '.fits', dustmap_11, dtype=['float32','float32','float32'], overwrite=True)
# print("d11 write file finished", "\n")