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]:





(array([  0.00000000e+00,   0.00000000e+00,   8.45770186e-05, ...,
          4.98854362e+00,   5.00760950e+00,   5.02675206e+00]),
 array([  0.00000000e+00,   0.00000000e+00,   2.90905674e-04, ...,
          2.08171820e+01,   2.08967455e+01,   2.09766292e+01]))

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]:





[<matplotlib.lines.Line2D at 0x7f0ccf1392d0>]






    

In [ ]:

    




    

In [ ]: