import fgbuster.mixingmatrix as mm
import healpy as hp
import numpy as np
from fgbuster import Dust, Synchrotron
from qubic.data import PATH as data_dir
[docs]
class CMBModel:
"""
CMB description assuming parametrized emission law such as :
Dl_CMB = r * Dl_tensor_r1 + Alens * Dl_lensed
Parameters
-----------
- params : Dictionary coming from `params.yml` file that define every parameters
- ell : Multipole used during the analysis
"""
def __init__(self, ell):
self.ell = ell
[docs]
def give_cl_cmb(self, r, Alens):
"""
Method to get theoretical CMB BB power spectrum according to Alens and r.
"""
power_spectrum = hp.read_cl(data_dir + "Cls_Planck2018_lensed_scalar.fits")[:, :4000]
if Alens != 1.0:
power_spectrum *= Alens
if r:
power_spectrum += r * hp.read_cl(data_dir + "Cls_Planck2018_unlensed_scalar_and_tensor_r1.fits")[:, :4000]
return power_spectrum
[docs]
def cl2dl(self, ell, cl):
"""
Method to convert Cl to Dl which is Dl = ell * (ell + 1) * Cl / 2 * pi
Arguments :
-----------
- ell : Array containing multipoles.
- cl : Array containing BB power spectrum.
"""
dl = np.zeros(ell.shape[0])
for i in range(ell.shape[0]):
dl[i] = (ell[i] * (ell[i] + 1) * cl[i]) / (2 * np.pi)
return dl
[docs]
def get_Dl_cmb(self):
"""
Method to interpolate theoretical BB power spectrum for effective multipoles.
"""
allDl = self.cl2dl(np.arange(1, 4001, 1), self.give_cl_cmb()[2])
Dl_eff = np.interp(self.ell, np.arange(1, 4001, 1), allDl)
return Dl_eff
[docs]
class SkySpectra:
def __init__(self, ell, nus, nu0_d=353, nu0_s=23, comp=None):
self.nus = nus
self.ell = ell
if self.ell is None:
self.ell = np.arange(2, 4000)
self.nbins = len(self.ell)
self.nspectra = len(self.nus)
self.nu0_d = nu0_d
self.nu0_s = nu0_s
if comp is not None:
self.comp = comp
self.nspectra = len(self.comp)
[docs]
def cl_to_dl(self, cl):
"""
Function to convert the cls into the dls
"""
return (self.ell * (self.ell + 1) * cl) / (2 * np.pi)
def _get_cl_cmb(self, r, Alens):
"""
Function to compute the CMB power spectrum from the Planck data
"""
power_spectrum = hp.read_cl(data_dir + "Cls_Planck2018_lensed_scalar.fits")[:, :4000]
if Alens != 1.0:
power_spectrum[2] *= Alens
if r:
power_spectrum += r * hp.read_cl(data_dir + "Cls_Planck2018_unlensed_scalar_and_tensor_r1.fits")[:, :4000]
return np.interp(self.ell, np.linspace(1, 4001, 4000), power_spectrum[2])
[docs]
def model_cmb(self, r, Alens):
"""
Define the CMB model, depending on r and Alens
"""
models = np.zeros((self.nspectra, self.nspectra, len(self.ell)))
cmb_ps = self._get_cl_cmb(r, Alens)
models = np.tile(cmb_ps, (self.nspectra, self.nspectra, 1))
return self.cl_to_dl(models)
[docs]
def scale_dust(self, betad, temp=20):
"""
Function to compute the dust mixing matrix element, depending on the frequency
"""
comp = Dust(nu0=self.nu0_d, temp=temp, beta_d=betad)
A = mm.MixingMatrix(comp).eval(self.nus)[:, 0]
return A[None, :] * A[:, None]
[docs]
def scale_sync(self, betas):
"""
Function to compute the dust mixing matrix element, depending on the frequency
"""
comp = Synchrotron(nu0=self.nu0_s, beta_pl=betas)
A = mm.MixingMatrix(comp).eval(self.nus)[:, 0]
return A[None, :] * A[:, None]
[docs]
def scale_dustsync(self, betad, betas, temp=20):
comp = Dust(nu0=self.nu0_d, temp=temp, beta_d=betad)
Adust = mm.MixingMatrix(comp).eval(self.nus)[:, 0]
comp = Synchrotron(nu0=self.nu0_s, beta_pl=betas)
Async = mm.MixingMatrix(comp).eval(self.nus)[:, 0]
return Adust[None, :] * Async[:, None] + Adust[:, None] * Async[None, :]
[docs]
def model(self, r, Alens, Ad, alphad, betad, As, alphas, betas, eps):
Dl_model = np.zeros((self.nspectra, self.nspectra, self.nbins))
### CMB
Dl_model += self.model_cmb(r, Alens)
### Foregrounds
prod_Anu_d = self.scale_dust(betad=betad)[..., None]
prod_Anu_s = self.scale_sync(betas=betas)[..., None]
Dl_model += Ad * prod_Anu_d * (self.ell / 80) ** alphad + As * prod_Anu_s * (self.ell / 80) ** alphas
prod_Anu_ds = self.scale_dustsync(betad=betad, betas=betas)[..., None]
Dl_model += eps * np.sqrt(abs(As) * abs(Ad)) * prod_Anu_ds * (self.ell / 80) ** ((alphad + alphas) / 2)
return Dl_model
[docs]
def modelCMM(self, r, Alens, Ad, alphad, betad, As, alphas, betas, eps):
# TODO : remove unnused args
Dl_model = np.zeros((self.nspectra, self.nspectra, self.nbins))
for i in range(self.nspectra):
if self.comp[i] == "CMB":
Dl_model[i, i] = self.model_cmb(r, Alens)[0, 0]
elif self.comp[i] == "Dust":
Dl_model[i, i] = Ad * (self.ell / 80) ** alphad
elif self.comp[i] == "Sync":
Dl_model[i, i] = As * (self.ell / 80) ** alphas
else:
raise ValueError(f"The component model {self.comp[i]} is not implement. Please use a known component.")
return Dl_model