In [2]:
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
In [9]:
from classy import Class
import numpy as np
import matplotlib.pyplot as plt
In [49]:
l_min = 2
l_max = 5000
Tcmb = 2.726
parameters = []
values = []
cosmo = Class()
def get_cls( values=[] ):
# Define your cosmology (what is not specified will be set to CLASS default parameters)
params = {
'output': 'tCl,pCl,lCl',
'l_max_scalars': l_max,
'lensing': 'yes'}
# set up the dictionary
for param, mean_value in zip(parameters, values):
params[param] = mean_value
# do Cl computation with classy
cosmo.set(params)
cosmo.compute()
cls = cosmo.lensed_cl(l_max)
cosmo.struct_cleanup()
cosmo.empty()
# convert dimensionless C_l's to C_l in muK**2, like MontePython
unit_factor = (Tcmb*1.0e6)**2 # used
return cls['ell'], cls['tt']*unit_factor,\
cls['te']*unit_factor, cls['ee']*unit_factor
def compute_error( theta_fwhm, sigma_T, sigma_P, f_sky,\
num_channels ):
# adapted from Monte Python
# convert arcmin to radians
theta_fwhm = theta_fwhm * np.array([np.pi/60./180.])
sigma_T = sigma_T * np.array([np.pi/60./180.])
sigma_P = sigma_P * np.array([np.pi/60./180.])
num_channels = num_channels
f_sky = f_sky
# compute noise in muK**2
noise_T = np.zeros(l_max+1, 'float64')
noise_P = np.zeros(l_max+1, 'float64')
for l in range(l_min, l_max+1):
noise_T[l] = 0
noise_P[l] = 0
for channel in range(num_channels):
noise_T[l] += sigma_T[channel]**-2 *\
np.exp(
-l*(l+1)*theta_fwhm[channel]**2/8./np.log(2.))
noise_P[l] += sigma_P[channel]**-2 *\
np.exp(
-l*(l+1)*theta_fwhm[channel]**2/8./np.log(2.))
noise_T[l] = 1/noise_T[l]
noise_P[l] = 1/noise_P[l]
return noise_T, noise_P
def get_DeltaCl( l, cl, nl, f_sky ):
return np.sqrt(2./((2*l+1)*f_sky)) * (cl + nl)
In [50]:
compute_error( num_channels = 3, f_sky = 0.65, \
theta_fwhm = [10., 7.1, 5.0], \
sigma_T = [68.1, 42.6, 65.4],
sigma_P = [109.4, 81.3, 133.6] )
Out[50]:
In [ ]:
In [89]:
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12,8))
# plt.figure()
Nl_T, Nl_P = compute_error( num_channels = 3, f_sky = 0.65, \
theta_fwhm = [10., 7.1, 5.0], \
sigma_T = [68.1, 42.6, 65.4],
sigma_P = [109.4, 81.3, 133.6])
ell, Cl_TT, Cl_TE, Cl_EE = get_cls()
dCl = get_DeltaCl( ell, Cl_TT, np.ones(len(ell))*1e-4, 0.65 )
Dl_fac = ell*(ell+1)/2./np.pi
ax1.plot( ell[ell>l_min], Dl_fac[ell>l_min] * Cl_TT[ell>l_min], '-', label='$C_l^{TT}$' )
# plt.plot( ell[ell>l_min], Dl_fac[ell>l_min] * dCl[ell>l_min], '-',
# label=r'$\Delta C_l$ with $N_{\ell} = 10^{-4} \, \mu \text{K}/T^2$' )
dCl_Planck = get_DeltaCl( ell, Cl_TT, Nl_T, 0.65 )
ax1.plot( ell[ell>l_min], Dl_fac[ell>l_min] * dCl_Planck[ell>l_min], '-',
label=r'Planck $\Delta C_l$ with white noise $N_{\ell}$' )
Dl_fac = ell*(ell+1)/2./np.pi
# plt.plot( ell[ell>l_min], Dl_fac[ell>l_min] * Cl_TT[ell>l_min], '-', label='$C_l^{TT}$' )
dCl_S4 = get_DeltaCl( ell, Cl_TT, np.ones(len(ell))*1e-4, 0.4 )
# plt.plot( ell[ell>l_min], Dl_fac[ell>l_min] * dCl_S4[ell>l_min], '-',
# label=r'$\Delta C_l$ with $N_{\ell} = 10^{-4} \, \mu \text{K}/T^2$' )
Nl_T, Nl_P = compute_error( num_channels = 1, f_sky = 0.4, \
theta_fwhm = [3.0], \
sigma_T = [1.0],
sigma_P = [1.4])
dCl_S4 = get_DeltaCl( ell, Cl_TT, Nl_T, 0.4 )
ax1.plot( ell[ell>l_min], Dl_fac[ell>l_min] * dCl_S4[ell>l_min], '-',
label=r'S4 $\Delta C_l$ with white noise $N_{\ell}$' )
ax1.set_yscale('log')
ax1.set_xscale('log')
ax1.set_xlim(1,3e3)
ax1.legend()
ax1.set_ylabel(r'$(l(l+1)/2\pi) C_l^{TT}$')
ax1.set_xlabel('l')
ax1.set_title('Planck $C_l^{TT}$ vs. Errors')
ax2.set_yscale('log')
ax2.set_xscale('log')
# ax1.set_xlim(1,3e3)
# ax2.legend()
ax2.set_ylabel(r'$(l(l+1)/2\pi) C_l^{TT}$')
ax2.set_xlabel('l')
ax2.set_title(r'$\Delta C_l(\text{S4}) / \Delta C_l(\text{Planck})$')
ax2.plot( ell[ell>l_min], (dCl_S4[ell>l_min] / dCl_Planck[ell>l_min]), '-',
label=r'S4 $\Delta C_l$ with white noise $N_{\ell}$' )
Out[89]:
In [ ]:
In [ ]: