Initialize

# import necessary modules
# uncomment to get plots displayed in notebook
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from classy_sz import Class as Class_sz
import math

font = {'size'   : 16, 'family':'STIXGeneral'}
axislabelfontsize='large'
matplotlib.rc('font', **font)
plt.rcParams["figure.figsize"] = [8.0,8.0]


plt.rcParams.update({
    "text.usetex": True,
    "font.family": "sans-serif",
    "font.sans-serif": ["Helvetica"]})

Maniyar et al Model

Compute

# a simple conversion from cl's to dl's
def cl_to_dl(lp):
    return lp*(lp+1.)/2./np.pi


Omegam0 = 0.3075
H0 = 67.74
Omegab = 0.0486
Omegac = 0.2589
Omegac + Omegab
omega_b = Omegab*(H0/100.)**2
omega_c = Omegac*(H0/100.)**2
hparam = H0/100.
maniyar_cosmo = {
'omega_b': omega_b,
'omega_cdm':  omega_c,
'h': H0/100.,
'ln10^{10}A_s': 3.048,
'n_s': 0.9665,

'm_ncdm': 0.0,
'cosmo_model': 1, # set to 1 for mnu-LCDM emulators and set mnu to 0.
}

%%time
class_sz = Class_sz()
class_sz.set({'output': 'lens_cib_1h,lens_cib_2h'})
class_sz.set(maniyar_cosmo)
class_sz.set({

'mass_function' : 'T08M200c',
'use_maniyar_cib_model':1,

'maniyar_cib_etamax' : 5.12572945e-01,

'maniyar_cib_zc' : 1.5,
'maniyar_cib_tau' : 8.25475287e-01,
'maniyar_cib_fsub' : 0.134*np.log(10.),
'Most_efficient_halo_mass_in_Msun' : 5.34372069e+12,
'Size_of_halo_masses_sourcing_CIB_emission' :  1.5583436676980493,
#for the Lsat tabulation:
'freq_min': 9e1,
'freq_max': 8.57e2,
'dlogfreq' : 0.1,

'concentration_parameter':'fixed', # this sets it to 5

'n_z_L_sat' :100,
'n_m_L_sat' :100,
'n_nu_L_sat':100,

'use_nc_1_for_all_halos_cib_HOD': 1,

'sub_halo_mass_function' : 'TW10',#'JvdB14',
'M_min_subhalo_in_Msun' : 1e5, # 1e5 see https://github.com/abhimaniyar/halomodel_cib_tsz_cibxtsz/blob/master/Cell_cib.py
'use_redshift_dependent_M_min': 0,
'M_min' : 1e8*hparam,
'M_max' : 1e15*hparam,
'z_min' : 0.012,
'z_max' : 10.,
'ell_min': 10.,
'ell_max':5e4,
'dlogell':0.3,


'ndim_redshifts': 210,
'ndim_masses':150,

'has_cib_flux_cut': 0,
'hm_consistency':0,

'epsabs_L_sat': 1e-40,
'epsrel_L_sat': 1e-9,

'damping_1h_term':0,

})

class_sz.set({
       'cib_frequency_list_num' : 6,
       'cib_frequency_list_in_GHz' : '100,143,217,353,545,857',
      })

CPU times: user 105 µs, sys: 155 µs, total: 260 µs
Wall time: 701 µs

True

%%time
class_sz.compute_class_szfast()

/Users/boris/Work/CLASS-SZ/SO-SZ/mcfit/mcfit/mcfit.py:130: UserWarning: use backend='jax' if desired
  warnings.warn("use backend='jax' if desired")

CPU times: user 2min 22s, sys: 4.81 s, total: 2min 27s
Wall time: 22.4 s

Dl_lens_cib = class_sz.cl_lens_cib()

Plot

nu_list = class_sz.pars['cib_frequency_list_in_GHz'].split(',')
# nu_list

from matplotlib.lines import Line2D
Planck = {'name': 'Planck',
          'do_cib': 1, 'do_tsz': 1, 'do_cibxtsz': 1,
          'freq_cib': [100., 143., 217., 353., 545., 857.],
          'cc': np.array([1.076, 1.017, 1.119, 1.097, 1.068, 0.995, 0.960]),
          'cc_cibmean': np.array([1.076, 1.017, 1.119, 1.097, 1.068, 0.995, 0.960]),
          'freq_cibmean': np.array([100., 143., 217., 353., 545., 857.]),
          'fc': np.ones(7),
          }
#Get Data
ells = np.asarray(Dl_lens_cib['545']['ell'])
factor = cl_to_dl(ells)

#Font Settings
label_size = 15
title_size = 18
legend_size = 13


fig, (ax1) = plt.subplots(1,1,figsize=(10,5))
ax = ax1
ax.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)
ax.tick_params(axis = 'y',which='both',length=5,direction='in', pad=5)
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
plt.setp(ax.get_yticklabels(), rotation='horizontal', fontsize=label_size)
plt.setp(ax.get_xticklabels(), fontsize=label_size)

#Colors
colors = ['k','r','b','orange','green','magenta']
alphas = [1,1,1,1,1,1]

for i, nu in enumerate(nu_list):
    #Extract Data
    nu_name = str(nu)
    cl_1h = np.asarray(Dl_lens_cib[str(nu)]['1h']) / factor
    cl_2h = np.asarray(Dl_lens_cib[str(nu)]['2h']) / factor
    cl_tot = cl_1h + cl_2h

    #Plot
    plt.plot(ells, ells * cl_tot, label= nu_name + ' GHz', color=colors[i],alpha= alphas[i])
    plt.plot(ells, ells * cl_1h, '-.', alpha= alphas[i], color=colors[i],lw=0.8)
    plt.plot(ells, ells * cl_2h, '--', alpha= alphas[i], color=colors[i],lw=0.8)

    #Figure Size
    # ax = plt.gca()
    # fig = plt.gcf()
    # fig.set_figheight(6)
    # fig.set_figwidth(11)


    #Legend
    line_1h = Line2D([0], [0], label=r'$\mathrm{1h}$', color='black', linestyle='-.', lw=0.8)
    line_2h = Line2D([0], [0], label=r'$\mathrm{2h}$', color='black', linestyle='--',lw=0.8)
    line_1p2h = Line2D([0], [0], label=r'$\mathrm{Total}$', color='black', linestyle='-')
    h, l = ax.get_legend_handles_labels()
    h.extend([line_1h, line_2h,line_1p2h])
    plt.legend(handles= h,frameon=False,fontsize=15,loc=2)

    #Labels
    plt.xlabel(r'$\ell$', fontsize= title_size,labelpad=1)
    plt.ylabel(r'$\ell C^{\nu\kappa}_{\ell}$\quad$\mathrm{[Jy/sr]}$', size= title_size)
    plt.xscale('log')
    plt.yscale('log')

    #Ticks
    # ax.yaxis.set_minor_locator(matpotlib.ticker.MultipleLocator(0.5))
    ax.tick_params(axis = 'x',which='both',length=5,direction='in', top=True, right=True, labelsize=label_size, pad=10)
    ax.tick_params(axis = 'y',which='both',length=5,direction='in', top=True, right=True, labelsize=label_size, pad=5)
    ax.grid( visible=True, which="both", alpha=0.2, linestyle='--')

plt.tight_layout()
# plt.savefig('figures/cls_CIB_CMB_lensing.pdf')

McCarthy & Madavacheril Model

import numpy as np
import matplotlib as mpl
import matplotlib
from matplotlib import pyplot as plt
from matplotlib_inline.backend_inline import set_matplotlib_formats
set_matplotlib_formats('svg')
from classy_sz import Class as Class_sz
import os
import time

font = {'family':'STIXGeneral'}
axislabelfontsize='large'
matplotlib.rc('font', **font)

plt.rcParams.update({
    "text.usetex": True,
    "font.family": "sans-serif",
    "font.sans-serif": ["Helvetica"]})

# a simple conversion from cl's to dl's
def cl_to_dl(lp):
    return lp*(lp+1.)/2./np.pi

You can set class_sz parameters individually, in a single dictionary, or multiple dictionaries.

First, let’s establish the cosmology. McCarthy et al. use Planck 2014 cosmology (https://arxiv.org/pdf/1303.5076.pdf, best-fit).

#Parameters for Cosmology Planck 14, https://arxiv.org/pdf/1303.5076.pdf, best-fit
p14_dict={}
p14_dict['h'] = 0.6711
p14_dict['omega_b'] = 0.022068
p14_dict['Omega_cdm'] = 0.3175-0.022068/0.6711/0.6711
p14_dict['A_s'] = 2.2e-9
p14_dict['n_s'] = .9624
# p14_dict['k_pivot'] = 0.05
p14_dict['tau_reio'] = 0.0925
# p14_dict['N_ncdm'] = 1
# p14_dict['N_ur'] = 0.00641
# p14_dict['deg_ncdm'] = 3
# p14_dict['m_ncdm'] = 0.02
# p14_dict['T_ncdm'] = 0.71611
p14_dict['cosmo_model'] = 1

Halo Model and grid parameters:

p_hm_dict = {}

p_hm_dict['mass_function'] = 'T10M200m'
p_hm_dict['concentration_parameter'] = 'D08'
p_hm_dict['delta_for_cib'] = '200m'
p_hm_dict['hm_consistency'] = 1
p_hm_dict['damping_1h_term'] = 0


# HOD parameters for CIB
p_hm_dict['M_min_HOD'] = pow(10.,10)  # was M_min_HOD_cib

#Grid Parameters
# Mass bounds
p_hm_dict['M_min'] = 1e8 * p14_dict['h']    # was M_min_cib
p_hm_dict['M_max'] = 1e16 * p14_dict['h']   # was M_max_cib
# Redshift bounds
p_hm_dict['z_min'] = 0.07
p_hm_dict['z_max'] = 6. # fiducial for MM20 : 6
p_hm_dict['freq_min'] = 10.
p_hm_dict['freq_max'] = 5e4
p_hm_dict['z_max_pk'] = p_hm_dict['z_max']

#Precision Parameters
# Precision for redshift integal
p_hm_dict['redshift_epsabs'] = 1e-40#1.e-40
p_hm_dict['redshift_epsrel'] = 1e-4#1.e-10 # fiducial value 1e-8
# Precision for mass integal
p_hm_dict['mass_epsabs'] = 1e-40 #1.e-40
p_hm_dict['mass_epsrel'] = 1e-4#1e-10
# Precision for Luminosity integral (sub-halo mass function)
p_hm_dict['L_sat_epsabs'] = 1e-40 #1.e-40
p_hm_dict['L_sat_epsrel'] = 1e-3#1e-10
# Multipole array
p_hm_dict['dlogell'] = 0.1
p_hm_dict['ell_max'] = 8000.
p_hm_dict['ell_min'] = 10

CIB Model parameters (see McCarthy et al. at https://journals.aps.org/prd/abstract/10.1103/PhysRevD.103.103515):

#CIB Parameters
p_CIB_dict = {}
p_CIB_dict['Redshift_evolution_of_L_-_M_normalisation'] =  0.36
p_CIB_dict['Dust_temperature_today_in_Kelvins'] = 24.4
p_CIB_dict['Emissivity_index_of_sed'] = 1.75
p_CIB_dict['Power_law_index_of_SED_at_high_frequency'] = 1.7
p_CIB_dict['Redshift_evolution_of_L_-_M_normalisation'] = 3.6
p_CIB_dict['Most_efficient_halo_mass_in_Msun'] = 10**12.6
p_CIB_dict['Normalisation_of_L_-_M_relation_in_[JyMPc2/Msun]'] = 6.4e-8
p_CIB_dict['Size_of_halo_masses_sourcing_CIB_emission'] = 0.5

Select your frequencies. Here, we choose some Planck frequencies. Note that the frequency lists (both for the observed CIB frequencies and the flux cuts) must be a single string, not a list of strings.

nu_list = [353,545,857]
nu_list_str = str(nu_list)[1:-1]  # Note: this must be a single string, not a list of strings!

#Frequency Parameters
p_freq_dict = {}
p_freq_dict['cib_frequency_list_num'] = len(nu_list)
p_freq_dict['cib_frequency_list_in_GHz'] = nu_list_str

We choose flux cuts from Planck 2013 (see https://journals.aps.org/prd/abstract/10.1103/PhysRevD.103.103515).

#Flux Cuts
cib_fcut_dict = {}

#Planck flux cut, Table 1 in https://arxiv.org/pdf/1309.0382.pdf
cib_fcut_dict['100'] = 400
cib_fcut_dict['143'] = 350
cib_fcut_dict['217'] = 225
cib_fcut_dict['353'] = 315
cib_fcut_dict['545'] = 350
cib_fcut_dict['857'] = 710
cib_fcut_dict['3000'] = 1000

def make_flux_cut_list(cib_flux, nu_list):
    """
    Make a string of flux cut values for given frequency list to pass into class_sz
    Beware: if frequency not in the flux_cut dictionary, it assigns 0
    """
    cib_flux_list = []
    keys = list(cib_flux.keys())
    for i,nu in enumerate(nu_list):
        if str(nu) in keys:
            cib_flux_list.append(cib_flux[str(nu)])
        else:
            cib_flux_list.append(0)
    return cib_flux_list

cib_flux_list = make_flux_cut_list(cib_fcut_dict, nu_list)

#Add Flux Cuts
p_freq_dict['cib_Snu_cutoff_list_in_mJy'] = str(list(cib_flux_list))[1:-1]
p_freq_dict['has_cib_flux_cut'] = 1

Compute

class_sz = Class_sz()

#Set Parameters
class_sz.set(p14_dict)
class_sz.set(p_hm_dict)
class_sz.set(p_CIB_dict)
class_sz.set(p_freq_dict)

#Tell It What You Want
class_sz.set({
    'output': 'lens_cib_1h,lens_cib_2h'
})

True

%%time
#Calculate Spectra
class_sz.set({
# 'ndim_redshifts':150,
})
class_sz.compute_class_szfast()
Dl_lens_cib = class_sz.cl_lens_cib()

/Users/boris/Work/CLASS-SZ/SO-SZ/mcfit/mcfit/mcfit.py:130: UserWarning: use backend='jax' if desired
  warnings.warn("use backend='jax' if desired")

CPU times: user 36.8 s, sys: 4.53 s, total: 41.3 s
Wall time: 7.75 s

Plot

#Get Data
ells = np.asarray(Dl_lens_cib['545']['ell'])
factor = cl_to_dl(ells)

#Font Settings
label_size = 15
title_size = 18
legend_size = 13


fig, (ax1) = plt.subplots(1,1,figsize=(10,5))
ax = ax1
ax.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)
ax.tick_params(axis = 'y',which='both',length=5,direction='in', pad=5)
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
plt.setp(ax.get_yticklabels(), rotation='horizontal', fontsize=label_size)
plt.setp(ax.get_xticklabels(), fontsize=label_size)

#Colors
colors = ['k','r','b']
alphas = [1,0.6,0.4]

for i, nu in enumerate(nu_list):
    #Extract Data
    nu_name = str(nu)
    cl_1h = np.asarray(Dl_lens_cib[str(nu)]['1h']) / factor
    cl_2h = np.asarray(Dl_lens_cib[str(nu)]['2h']) / factor
    cl_tot = cl_1h + cl_2h

    #Plot
    plt.plot(ells, ells * cl_tot, label= nu_name + ' GHz', color=colors[i],alpha= alphas[i])
    plt.plot(ells, ells * cl_1h, '-.', alpha= alphas[i], color=colors[i],lw=0.8)
    plt.plot(ells, ells * cl_2h, '--', alpha= alphas[i], color=colors[i],lw=0.8)

    #Figure Size
    # ax = plt.gca()
    # fig = plt.gcf()
    # fig.set_figheight(6)
    # fig.set_figwidth(11)


    #Legend
    line_1h = mpl.lines.Line2D([0], [0], label=r'$\mathrm{1h}$', color='black', linestyle='-.', lw=0.8)
    line_2h = mpl.lines.Line2D([0], [0], label=r'$\mathrm{2h}$', color='black', linestyle='--',lw=0.8)
    line_1p2h = mpl.lines.Line2D([0], [0], label=r'$\mathrm{Total}$', color='black', linestyle='-')
    h, l = ax.get_legend_handles_labels()
    h.extend([line_1h, line_2h,line_1p2h])
    plt.legend(handles= h,frameon=False,fontsize=15,loc=2)

    #Labels
    plt.xlabel(r'$\ell$', fontsize= title_size,labelpad=1)
    plt.ylabel(r'$\ell C^{\nu\kappa}_{\ell}$\quad$\mathrm{[Jy/sr]}$', size= title_size)
    plt.xscale('log')
    # plt.yscale('log')

    #Ticks
    ax.yaxis.set_minor_locator(mpl.ticker.MultipleLocator(0.5))
    ax.tick_params(axis = 'x',which='both',length=5,direction='in', top=True, right=True, labelsize=label_size, pad=10)
    ax.tick_params(axis = 'y',which='both',length=5,direction='in', top=True, right=True, labelsize=label_size, pad=5)
    ax.grid( visible=True, which="both", alpha=0.2, linestyle='--')

plt.tight_layout()
# plt.savefig('figures/cls_CIB_CMB_lensing.pdf')

Primordial non-Gaussianity

Initialize cosmology part

%%time
class_sz = Class_sz()
#Set Parameters
class_sz.set(p14_dict)
class_sz.set(p_hm_dict)
class_sz.set(p_CIB_dict)
class_sz.set(p_freq_dict)
#Calculate Spectra
class_sz.set({
# 'ndim_redshifts':150,
'use_png_scale_dependent_bias' : 1,
'fNL' : 100.,
'hm_consistency': 1,
'ndim_redshifts': 100,
'ndim_masses': 100,
})
#Tell It What You Want
class_sz.set({
    'output': 'lens_cib_1h,lens_cib_2h'
})
class_sz.compute()

CPU times: user 1min 2s, sys: 600 ms, total: 1min 3s
Wall time: 10.5 s

Compute (without redoing cosmo) but only changing fNL:

%%time
## define array of fNL
fnl_ar = np.linspace(-100,100,10)
Dl_lens_cib_all = []
for fnl in fnl_ar:
    class_sz.compute_class_sz({'fNL':fnl})
    Dl_lens_cib_all.append(class_sz.cl_lens_cib())

CPU times: user 4min 35s, sys: 1.17 s, total: 4min 36s
Wall time: 40 s

#Get Data
ells = np.asarray(Dl_lens_cib['545']['ell'])
factor = cl_to_dl(ells)

#Font Settings
label_size = 15
title_size = 18
legend_size = 13


fig, (ax1) = plt.subplots(1,1,figsize=(10,5))
ax = ax1
ax.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)
ax.tick_params(axis = 'y',which='both',length=5,direction='in', pad=5)
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
plt.setp(ax.get_yticklabels(), rotation='horizontal', fontsize=label_size)
plt.setp(ax.get_xticklabels(), fontsize=label_size)

#Colors
colors = ['k','r','b']
alphas = [1,0.6,0.4]
for ifnl,fnl in enumerate(fnl_ar):
    Dl_lens_cib = Dl_lens_cib_all[ifnl]
    for i, nu in enumerate(nu_list):
        #Extract Data
        nu_name = str(nu)
        cl_1h = np.asarray(Dl_lens_cib[str(nu)]['1h']) / factor
        cl_2h = np.asarray(Dl_lens_cib[str(nu)]['2h']) / factor
        cl_tot = cl_1h + cl_2h

        #Plot
        plt.plot(ells, ells * cl_tot, label= nu_name + ' GHz', color=colors[i],alpha= alphas[i])
        plt.plot(ells, ells * cl_1h, '-.', alpha= alphas[i], color=colors[i],lw=0.8)
        plt.plot(ells, ells * cl_2h, '--', alpha= alphas[i], color=colors[i],lw=0.8)

        #Figure Size
        # ax = plt.gca()
        # fig = plt.gcf()
        # fig.set_figheight(6)
        # fig.set_figwidth(11)


        #Legend
        line_1h = mpl.lines.Line2D([0], [0], #label=r'$\mathrm{1h}$',
                                   color='black', linestyle='-.', lw=0.8)
        line_2h = mpl.lines.Line2D([0], [0], #label=r'$\mathrm{2h}$',
                                   color='black', linestyle='--',lw=0.8)
        line_1p2h = mpl.lines.Line2D([0], [0], #label=r'$\mathrm{Total}$',
                                     color='black', linestyle='-')
        # h, l = ax.get_legend_handles_labels()
        h.extend([line_1h, line_2h,line_1p2h])
        # plt.legend(handles= h,frameon=False,fontsize=15,loc=2)

        #Labels
        plt.xlabel(r'$\ell$', fontsize= title_size,labelpad=1)
        plt.ylabel(r'$\ell C^{\nu\kappa}_{\ell}$\quad$\mathrm{[Jy/sr]}$', size= title_size)
        plt.xscale('log')
        # plt.yscale('log')

        #Ticks
        ax.yaxis.set_minor_locator(mpl.ticker.MultipleLocator(0.5))
        ax.tick_params(axis = 'x',which='both',length=5,direction='in', top=True, right=True, labelsize=label_size, pad=10)
        ax.tick_params(axis = 'y',which='both',length=5,direction='in', top=True, right=True, labelsize=label_size, pad=5)
        ax.grid( visible=True, which="both", alpha=0.2, linestyle='--')

plt.tight_layout()
# plt.savefig('figures/cls_CIB_CMB_lensing.pdf')