In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import sys, platform, os
import matplotlib
from matplotlib import pyplot as plt
import numpy as np

#Assume installed from github using "git clone --recursive https://github.com/cmbant/CAMB.git"
#This file is then in the docs folders
isitgr_path = os.path.realpath(os.path.join(os.getcwd(),'..'))
sys.path.insert(0,isitgr_path)
import isitgr
from isitgr import model, initialpower
print('Using CAMB-ISiTGR %s installed at %s'%(isitgr.__version__,os.path.dirname(isitgr.__file__)))

Using CAMB-ISiTGR 1.1.0 installed at /home/cristhian/Projects/ISiTGR-2020/ISiTGR-2020/camb/isitgr

($\mu$, $\eta$) parameterization: Reproducing Fig. 1 of arXiv:1502.01590v2

(slight differences are present because cosmological parameter values are now a bit different)

In [2]:
# Defining array for different MG parameter values.
E11 = [1, -1, 0.5, 0]
E22 = [1, -1, 0.5, 1]

# Set up different set of parameters for CAMB (including MG)
pars_GR = isitgr.CAMBparams()
pars_MG1 = isitgr.CAMBparams()
pars_MG2 = isitgr.CAMBparams()
pars_MG3 = isitgr.CAMBparams()
pars_MG4 = isitgr.CAMBparams()

# This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency

# For GR
pars_GR.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09)
pars_GR.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_GR.set_for_lmax(2500, lens_potential_accuracy=0);

# For MG model 1. Here we set parameterization="mueta". 
# No scale dependence is considered in any MG model since we do not call c1, c2 or Lambda parameters...
pars_MG1.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="mueta", E11=E11[0], E22=E22[0])
pars_MG1.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG1.set_for_lmax(2500, lens_potential_accuracy=0);

# For MG model 2. Here we set parameterization="mueta". 
pars_MG2.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="mueta", E11=E11[1], E22=E22[1])
pars_MG2.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG2.set_for_lmax(2500, lens_potential_accuracy=0);

# For MG model 3. Here we set parameterization="mueta". 
pars_MG3.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="mueta", E11=E11[2], E22=E22[2])
pars_MG3.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG3.set_for_lmax(2500, lens_potential_accuracy=0);

# For MG model 4. Here we set parameterization="mueta". 
pars_MG4.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="mueta", E11=E11[3], E22=E22[3])
pars_MG4.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG4.set_for_lmax(2500, lens_potential_accuracy=0);
In [3]:
# Calculate results for different models.
results_GR = isitgr.get_results(pars_GR)
results_MG1 = isitgr.get_results(pars_MG1)
results_MG2 = isitgr.get_results(pars_MG2)
results_MG3 = isitgr.get_results(pars_MG3)
results_MG4 = isitgr.get_results(pars_MG4)

# (1) camb.get_results() compute results for specified parameters and return CAMBdata.
# (2) camb.results.CAMBdata is an object for storing calculational data, parameters and transfer functions.
In [4]:
# Get dictionary of CAMB power spectra for different models, including GR and MG models.
powers_GR =results_GR.get_cmb_power_spectra(pars_GR, CMB_unit='muK')
powers_MG1 =results_MG1.get_cmb_power_spectra(pars_MG1, CMB_unit='muK')
powers_MG2 =results_MG2.get_cmb_power_spectra(pars_MG2, CMB_unit='muK')
powers_MG3 =results_MG3.get_cmb_power_spectra(pars_MG3, CMB_unit='muK')
powers_MG4 =results_MG4.get_cmb_power_spectra(pars_MG4, CMB_unit='muK')

# For name in powers_GR: print(name)
In [5]:
# Here we obtain the total lensed CMB power spectra
totCL_GR=powers_GR['total']
totCL_MG1=powers_MG1['total']
totCL_MG2=powers_MG2['total']
totCL_MG3=powers_MG3['total']
totCL_MG4=powers_MG4['total']
print(totCL_GR.shape)

(2551, 4)
In [6]:
# Here we obtain the cl's for lensing potential.

# For GR
pars_GR.set_dark_energy(w=-1, wa=0, dark_energy_model='fluid') 
cl_GR = results_GR.get_lens_potential_cls(lmax=2550)

# For MG1
pars_MG1.set_dark_energy(w=-1, wa=0, dark_energy_model='fluid') 
cl_MG1 = results_MG1.get_lens_potential_cls(lmax=2550)

# For MG2
pars_MG2.set_dark_energy(w=-1, wa=0, dark_energy_model='fluid') 
cl_MG2 = results_MG2.get_lens_potential_cls(lmax=2550)

# For MG3
pars_MG3.set_dark_energy(w=-1, wa=0, dark_energy_model='fluid') 
cl_MG3 = results_MG3.get_lens_potential_cls(lmax=2550)

# For MG4
pars_MG4.set_dark_energy(w=-1, wa=0, dark_energy_model='fluid') 
cl_MG4 = results_MG4.get_lens_potential_cls(lmax=2550)
In [7]:
Tspectra = [totCL_GR[:,0], totCL_MG1[:,0], totCL_MG2[:,0], totCL_MG3[:,0], totCL_MG4[:,0]]
Lspectra = [cl_GR[:,0], cl_MG1[:,0], cl_MG2[:,0], cl_MG3[:,0], cl_MG4[:,0]]
In [8]:
from matplotlib import gridspec
import matplotlib.ticker as mtick
# Desired overall figure size
SIZE = (10, 10)
# GridSpect(nrows,ncolumns,ratio)
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1]) 
# The figure
fig = plt.figure(figsize = SIZE)

y_values = [Tspectra, Lspectra]
x_values = np.arange(2551)

# Creating the subplots
base = 0
ax = []
for y,h in zip(y_values[0],[1,2]):
    ax.append(fig.add_subplot(gs[base:h+base,:]))
    base += h
    ax[-1].semilogx(x_values,y_values[-1+h][0], color='k', label='$\Lambda CDM$')
    ax[-1].semilogx(x_values,y_values[-1+h][1], color='b', label='$E_{11}=1$, $E_{22}=1$')
    ax[-1].semilogx(x_values,y_values[-1+h][2], color='cyan', linestyle='dashdot', label='$E_{11}=-1$, $E_{22}=-1$')
    ax[-1].semilogx(x_values,y_values[-1+h][3], color='teal', linestyle='dotted', label='$E_{11}=0.5$, $E_{22}=0.5$')
    ax[-1].semilogx(x_values,y_values[-1+h][4], color='darkblue', linestyle='dashed', label='$E_{11}=0$, $E_{22}=1$')

# Plot spacing (this is just to concatenate them)
fig.subplots_adjust(hspace=0.001)

# Some final settings
ax[0].set_title('obtained using ISiTGR-python wrapper')
ax[0].set_ylabel('$\ell(\ell+1)C_{\ell}^{TT}/2\pi$ [$\mu K^2$]', fontsize=16)
ax[1].set_ylabel('$[\ell(\ell+1)]^2C_{\ell}^{\phi\phi}/2\pi$ [$\mu K^2$]', fontsize=16)
plt.xlabel('$\ell$', fontsize=16)
ax[0].set_xlim([2, 2100])
ax[1].set_xlim([2, 2100])
ax[0].set_ylim([0, 8000])
ax[1].set_ylim([0, 4E-7])
ax[1].yaxis.set_major_formatter(mtick.FormatStrFormatter('%.e'))
ax[1].legend(loc='upper right', fontsize=16)
yticks = ax[0].yaxis.get_major_ticks() 
yticks[0].label1.set_visible(False)

($\mu$, $\Sigma$) parameterization: Reproducing Fig. 4 of arXiv:1212.3339v2

In [2]:
# Defining array for different MG parameter values.
mu0 = [0.5, 0]
Sigma0 = [0, -0.5]

# Set up different set of parameters for CAMB (including MG)
pars_GR = isitgr.CAMBparams()
pars_MG1 = isitgr.CAMBparams()
pars_MG2 = isitgr.CAMBparams()

# This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency

# For GR
pars_GR.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09)
pars_GR.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_GR.set_for_lmax(500, lens_potential_accuracy=0);

# For MG model 1. Here we set parameterization="mueta". 
# No scale dependence is considered in any MG model since we do not call c1, c2 or Lambda parameters...
pars_MG1.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="muSigma", mu0=mu0[0], Sigma0=Sigma0[0])
pars_MG1.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG1.set_for_lmax(500, lens_potential_accuracy=0);

# For MG model 2. Here we set parameterization="mueta". 
pars_MG2.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="muSigma", mu0=mu0[1], Sigma0=Sigma0[1])
pars_MG2.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG2.set_for_lmax(500, lens_potential_accuracy=0);
pars_MG1
Out[2]:
class: <CAMBparams>
 WantCls = True
 WantTransfer = False
 WantScalars = True
 WantTensors = False
 WantVectors = False
 WantDerivedParameters = True
 Want_cl_2D_array = True
 Want_CMB = True
 Want_CMB_lensing = True
 DoLensing = True
 NonLinear = NonLinear_none
 Transfer: <TransferParams>
   high_precision = False
   accurate_massive_neutrinos = False
   kmax = 0.9
   k_per_logint = 0
   PK_num_redshifts = 1
   PK_redshifts = [0.0]
 want_zstar = False
 want_zdrag = False
 min_l = 2
 max_l = 650
 max_l_tensor = 600
 max_eta_k = 1625.0
 max_eta_k_tensor = 1200.0
 ombh2 = 0.0226
 omch2 = 0.112
 omk = 0.0
 omnuh2 = 0.0
 H0 = 70.0
 TCMB = 2.7255
 YHe = 0.24548291592414984
 num_nu_massless = 3.046
 num_nu_massive = 0
 nu_mass_eigenstates = 0
 share_delta_neff = False
 nu_mass_degeneracies = []
 nu_mass_fractions = []
 nu_mass_numbers = []
 InitPower: <InitialPowerLaw>
   tensor_parameterization = tensor_param_rpivot
   ns = 0.96
   nrun = 0.0
   nrunrun = 0.0
   nt = 0.0
   ntrun = -0.0
   r = 0.0
   pivot_scalar = 0.05
   pivot_tensor = 0.05
   As = 2.1e-09
   At = 1.0
 Recomb: <Recfast>
   min_a_evolve_Tm = 0.0011098779505118728
   RECFAST_fudge = 1.125
   RECFAST_fudge_He = 0.86
   RECFAST_Heswitch = 6
   RECFAST_Hswitch = True
   AGauss1 = -0.14
   AGauss2 = 0.079
   zGauss1 = 7.28
   zGauss2 = 6.73
   wGauss1 = 0.18
   wGauss2 = 0.33
 Reion: <TanhReionization>
   Reionization = True
   use_optical_depth = True
   redshift = 10.0
   optical_depth = 0.09
   delta_redshift = 0.5
   fraction = -1.0
   include_helium_fullreion = True
   helium_redshift = 3.5
   helium_delta_redshift = 0.4
   helium_redshiftstart = 5.5
   tau_solve_accuracy_boost = 1.0
   timestep_boost = 1.0
   max_redshift = 50.0
 DarkEnergy: <DarkEnergyFluid>
   w = -1.0
   wa = 0.0
   cs2 = 1.0
   use_tabulated_w = False
 NonLinearModel: <Halofit>
   Min_kh_nonlinear = 0.005
   halofit_version = mead
   HMCode_A_baryon = 3.13
   HMCode_eta_baryon = 0.603
 Accuracy: <AccuracyParams>
   AccuracyBoost = 1.0
   lSampleBoost = 1.0
   lAccuracyBoost = 1.0
   AccuratePolarization = True
   AccurateBB = False
   AccurateReionization = True
   TimeStepBoost = 1.0
   BackgroundTimeStepBoost = 1.0
   IntTolBoost = 1.0
   SourcekAccuracyBoost = 1.0
   IntkAccuracyBoost = 1.0
   TransferkBoost = 1.0
   NonFlatIntAccuracyBoost = 1.0
   BessIntBoost = 1.0
   LensingBoost = 1.0
   NonlinSourceBoost = 1.0
   BesselBoost = 1.0
   LimberBoost = 1.0
   SourceLimberBoost = 1.0
   KmaxBoost = 1.0
   neutrino_q_boost = 1.0
 SourceTerms: <SourceTermParams>
   limber_windows = True
   limber_phi_lmin = 100
   counts_density = True
   counts_redshift = True
   counts_lensing = False
   counts_velocity = True
   counts_radial = False
   counts_timedelay = True
   counts_ISW = True
   counts_potential = True
   counts_evolve = False
   line_phot_dipole = False
   line_phot_quadrupole = False
   line_basic = True
   line_distortions = True
   line_extra = False
   line_reionization = False
   use_21cm_mK = True
 z_outputs = []
 scalar_initial_condition = initial_adiabatic
 InitialConditionVector = []
 OutputNormalization = 1
 Alens = 1.0
 MassiveNuMethod = Nu_best
 DoLateRadTruncation = True
 Evolve_baryon_cs = False
 Evolve_delta_xe = False
 Evolve_delta_Ts = False
 Do21cm = False
 transfer_21cm_cl = False
 Log_lvalues = False
 use_cl_spline_template = True
 SourceWindows = []
 CustomSources: <CustomSources>
   num_custom_sources = 0
   c_source_func = None
   custom_source_ell_scales = []
 Q0 = 0.0
 Qinf = 0.0
 DR0 = 0.0
 DRinf = 0.0
 s = 0.0
 k_c = 0.0
 Q1 = 0.0
 Q2 = 0.0
 Q3 = 0.0
 Q4 = 0.0
 D1 = 0.0
 D2 = 0.0
 D3 = 0.0
 D4 = 0.0
 z_div = 1.0
 z_TGR = 1.0
 z_tw = 1.0
 k_tw = 1.0
 t_dep = False
 R_func = False
 true_bin = False
 ISiTGR_BIN = False
 E11 = 0.0
 E22 = 0.0
 mu0 = 0.5
 Sigma0 = 0.0
 c1 = 1.0
 c2 = 1.0
 Lambda = 0.0
 mu1 = 0.0
 mu2 = 0.0
 mu3 = 0.0
 mu4 = 0.0
 eta1 = 0.0
 eta2 = 0.0
 eta3 = 0.0
 eta4 = 0.0
 Sigma1 = 0.0
 Sigma2 = 0.0
 Sigma3 = 0.0
 Sigma4 = 0.0
 ISiTGR_mueta = False
 ISiTGR_muSigma = True
 ISiTGR_QDR = False
 ISiTGR_BIN_mueta = False
 ISiTGR_BIN_muSigma = False
 GR = 1
In [3]:
# Calculate results for different models.
results_GR = isitgr.get_results(pars_GR)
results_MG1 = isitgr.get_results(pars_MG1)
results_MG2 = isitgr.get_results(pars_MG2)

# (1) camb.get_results() compute results for specified parameters and return CAMBdata.
# (2) camb.results.CAMBdata is an object for storing calculational data, parameters and transfer functions.
In [12]:
#get dictionary of CAMB power spectra for different models, including GR and MG models.
powers_GR =results_GR.get_cmb_power_spectra(pars_GR, CMB_unit='muK')
powers_MG1 =results_MG1.get_cmb_power_spectra(pars_MG1, CMB_unit='muK')
powers_MG2 =results_MG2.get_cmb_power_spectra(pars_MG2, CMB_unit='muK')
In [13]:
# Here we obtain the total lensed CMB power spectra
totCL_GR=powers_GR['total']
totCL_MG1=powers_MG1['total']
totCL_MG2=powers_MG2['total']
In [14]:
Tspectra = [totCL_GR[:,0], totCL_MG1[:,0], totCL_MG2[:,0]]
TEspectra = [totCL_GR[:,3], totCL_MG1[:,3], totCL_MG2[:,3]]
ell = x_values = np.arange(551)
TEspectra_divided_by_ell = TEspectra/ell #we divided first element by 0 so we get nan, but it's irrelevant

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in true_divide
  after removing the cwd from sys.path.
In [15]:
from matplotlib import gridspec
import matplotlib.ticker as mtick
# desired overall figure size
SIZE = (10, 15)
# GridSpect(nrows,ncolumns,ratio)
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1]) 
# The figure
fig = plt.figure(figsize = SIZE)

y_values = [Tspectra, TEspectra_divided_by_ell]

# creating the subplots
base = 0
ax = []
#for y,h in zip(y_values[0],[1,2]):
ax.append(fig.add_subplot(gs[base:1+base,:]))
ax[-1].loglog(ell,y_values[0][0], color='k', label='$\Lambda CDM$')
ax[-1].loglog(ell,y_values[0][1], color='b', linestyle='dashed', label='$\mu_{0}=0.5$, $\Sigma_{0}=0$')
ax[-1].loglog(ell,y_values[0][2], color='pink', linestyle='dashdot', label='$\mu_{0}=0$, $\Sigma_{0}=-0.5$')
ax.append(fig.add_subplot(gs[(base+1):1+(base+1),:]))
ax[-1].semilogx(ell,y_values[1][0], color='k', label='$\Lambda CDM$')
ax[-1].semilogx(ell,y_values[1][1], color='b', linestyle='dashed', label='$\mu_{0}=0.5$, $\Sigma_{0}=0$')
ax[-1].semilogx(ell,y_values[1][2], color='pink', linestyle='dashdot', label='$\mu_{0}=0$, $\Sigma_{0}=-0.5$')
# plot spacing (this is just to concatenate them)
fig.subplots_adjust(hspace=0.001)

# some final settings
ax[0].set_title('obtained using ISiTGR-python wrapper')
ax[0].set_ylabel('$\ell(\ell+1)C_{\ell}^{TT}/2\pi$ [$\mu K^2$]', fontsize=16)
ax[1].set_ylabel('$(\ell+1)C_{\ell}^{TE}/2\pi$ [$\mu K^2$]', fontsize=16)
plt.xlabel('$\ell$', fontsize=16)
ax[0].set_xlim([2, 500])
ax[1].set_xlim([2, 500])
ax[0].set_ylim([100, 100000])
ax[1].set_ylim([-1.5, 3])
ax[0].legend(loc='upper left', fontsize=16)
yticks = ax[1].yaxis.get_major_ticks() 
yticks[-1].label1.set_visible(False)

($Q$, $R$) parameterization: Reproducing top plot of Fig. 1 in arXiv:1002.4197v4

In [16]:
import math
# Defining array for different MG parameter values.
Q0 = [0.5, 1.5, 2]

# Set up different set of parameters for CAMB (including MG)
pars_GR = isitgr.CAMBparams()
pars_MG1 = isitgr.CAMBparams()
pars_MG2 = isitgr.CAMBparams()
pars_MG3 = isitgr.CAMBparams()

# This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency

# For GR
pars_GR.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09)
pars_GR.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_GR.set_for_lmax(1000, lens_potential_accuracy=0);

# For MG model 1. Here we set parameterization="mueta". 
# No scale dependence is considered in any MG model since we do not call c1, c2 or Lambda parameters...
pars_MG1.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="QR", Q0=Q0[0], s=3, k_c=math.inf, Qinf=0, R0=1, Rinf=1)
pars_MG1.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG1.set_for_lmax(1000, lens_potential_accuracy=0);

# For MG model 2. Here we set parameterization="mueta". 
pars_MG2.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="QR", Q0=Q0[1], s=3, k_c=math.inf, Qinf=0, R0=1, Rinf=1)
pars_MG2.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG2.set_for_lmax(1000, lens_potential_accuracy=0);

# For MG model 3. Here we set parameterization="mueta". 
pars_MG3.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                  parameterization="QR", Q0=Q0[2], s=3, k_c=math.inf, Qinf=0, R0=1, Rinf=1)
pars_MG3.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
pars_MG3.set_for_lmax(1000, lens_potential_accuracy=0);
In [17]:
# Calculate results for different models.
results_GR = isitgr.get_results(pars_GR)
results_MG1 = isitgr.get_results(pars_MG1)
results_MG2 = isitgr.get_results(pars_MG2)
results_MG3 = isitgr.get_results(pars_MG3)


# (1) camb.get_results() compute results for specified parameters and return CAMBdata.
# (2) camb.results.CAMBdata is an object for storing calculational data, parameters and transfer functions.
In [18]:
#get dictionary of CAMB power spectra for different models, including GR and MG models.
powers_GR =results_GR.get_cmb_power_spectra(pars_GR, CMB_unit='muK')
powers_MG1 =results_MG1.get_cmb_power_spectra(pars_MG1, CMB_unit='muK')
powers_MG2 =results_MG2.get_cmb_power_spectra(pars_MG2, CMB_unit='muK')
powers_MG3 =results_MG3.get_cmb_power_spectra(pars_MG3, CMB_unit='muK')
In [19]:
# Here we obtain the total lensed CMB power spectra
totCL_GR=powers_GR['total']
totCL_MG1=powers_MG1['total']
totCL_MG2=powers_MG2['total']
totCL_MG3=powers_MG3['total']
In [20]:
# Plot settings. We plot the CMB temperature power spectra
ls = np.arange(totCL_GR.shape[0])
plt.figure(figsize=(10,10))
plt.semilogx(ls,totCL_GR[:,0], color='k', label='$\Lambda CDM$')
plt.semilogx(ls,totCL_MG1[:,0], color='red', linestyle='dotted', label='$Q_0=0.5$, $s=3$')
plt.semilogx(ls,totCL_MG2[:,0], color='green', linestyle='dashdot', label='$Q_0=1.5$, $s=3$')
plt.semilogx(ls,totCL_MG3[:,0], color='blue', linestyle='dashed', label='$Q_0=2$, $s=3$')
plt.xlim([2,1000])
plt.ylim([0,6500])
plt.legend(loc='upper left', fontsize=16)
plt.ylabel('$\ell(\ell+1)C_{\ell}^{TT}/2\pi$ [$\mu K^2$]', fontsize=16)
plt.xlabel('$\ell$', fontsize=16)
plt.title('obtained using ISiTGR-python wrapper')
plt.text(12, 300, "$Q_{\infty}=0$, $R_{0}=1$, $R_{\infty}=1$, $k_c=\infty$.", fontsize=16)
Out[20]:
Text(12, 300, '$Q_{\\infty}=0$, $R_{0}=1$, $R_{\\infty}=1$, $k_c=\\infty$.')

Binning Method using ($Q$, $D$) parameterization: Reproducing Fig.2 in arXiv:1109.4583v3

In [21]:
# Set up different set of parameters for CAMB (including MG)
pars_GR = isitgr.CAMBparams()
pars_binning_hybrid = isitgr.CAMBparams()
pars_binning_trad = isitgr.CAMBparams()

# This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency

# For GR
pars_GR.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09)
pars_GR.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)

# For MG model 1. Here we set parameterization="mueta". 
# No scale dependence is considered in any MG model since we do not call c1, c2 or Lambda parameters...
pars_binning_hybrid.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                                  parameterization="QD", binning="hybrid", z_div=1, z_TGR=2, z_tw=0.05,
                                  k_c=0.01, Q1=1.5, Q2=0.5, Q3=1.5, Q4=0.5,
                                  D1=1, D2=1, D3=1, D4=1)
pars_binning_hybrid.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)

# For MG model 2. Here we set parameterization="mueta". 
pars_binning_trad.set_cosmology(H0=70, ombh2=0.0226, omch2=0.112, mnu=0, omk=0, tau=0.09,
                                  parameterization="QD", binning="traditional", z_div=1, z_TGR=2, z_tw=0.05,
                                  k_c=0.01, Q1=1.5, Q2=0.5, Q3=1.5, Q4=0.5,
                                  D1=1, D2=1, D3=1, D4=1, k_tw=0.001)
pars_binning_trad.InitPower.set_params(As=2.1e-9, ns=0.96, r=0)
Out[21]:
class: <InitialPowerLaw>
 tensor_parameterization = tensor_param_rpivot
 ns = 0.96
 nrun = 0.0
 nrunrun = 0.0
 nt = 0.0
 ntrun = -0.0
 r = 0.0
 pivot_scalar = 0.05
 pivot_tensor = 0.05
 As = 2.1e-09
 At = 1.0
In [22]:
#Note non-linear corrections couples to smaller scales than you want
pars_GR.set_matter_power(redshifts=(0,1), kmax=2)
pars_binning_hybrid.set_matter_power(redshifts=(0,1), kmax=2)
pars_binning_trad.set_matter_power(redshifts=(0,1), kmax=2)

Note: redshifts have been re-sorted (earliest first)
Note: redshifts have been re-sorted (earliest first)
Note: redshifts have been re-sorted (earliest first)
Out[22]:
class: <CAMBparams>
 WantCls = True
 WantTransfer = True
 WantScalars = True
 WantTensors = False
 WantVectors = False
 WantDerivedParameters = True
 Want_cl_2D_array = True
 Want_CMB = True
 Want_CMB_lensing = True
 DoLensing = True
 NonLinear = NonLinear_none
 Transfer: <TransferParams>
   high_precision = True
   accurate_massive_neutrinos = False
   kmax = 2.0
   k_per_logint = 0
   PK_num_redshifts = 2
   PK_redshifts = [1.0, 0.0]
 want_zstar = False
 want_zdrag = False
 min_l = 2
 max_l = 2500
 max_l_tensor = 600
 max_eta_k = 5000.0
 max_eta_k_tensor = 1200.0
 ombh2 = 0.0226
 omch2 = 0.112
 omk = 0.0
 omnuh2 = 0.0
 H0 = 70.0
 TCMB = 2.7255
 YHe = 0.24548291592414984
 num_nu_massless = 3.046
 num_nu_massive = 0
 nu_mass_eigenstates = 0
 share_delta_neff = False
 nu_mass_degeneracies = []
 nu_mass_fractions = []
 nu_mass_numbers = []
 InitPower: <InitialPowerLaw>
   tensor_parameterization = tensor_param_rpivot
   ns = 0.96
   nrun = 0.0
   nrunrun = 0.0
   nt = 0.0
   ntrun = -0.0
   r = 0.0
   pivot_scalar = 0.05
   pivot_tensor = 0.05
   As = 2.1e-09
   At = 1.0
 Recomb: <Recfast>
   min_a_evolve_Tm = 0.0011098779505118728
   RECFAST_fudge = 1.125
   RECFAST_fudge_He = 0.86
   RECFAST_Heswitch = 6
   RECFAST_Hswitch = True
   AGauss1 = -0.14
   AGauss2 = 0.079
   zGauss1 = 7.28
   zGauss2 = 6.73
   wGauss1 = 0.18
   wGauss2 = 0.33
 Reion: <TanhReionization>
   Reionization = True
   use_optical_depth = True
   redshift = 10.0
   optical_depth = 0.09
   delta_redshift = 0.5
   fraction = -1.0
   include_helium_fullreion = True
   helium_redshift = 3.5
   helium_delta_redshift = 0.4
   helium_redshiftstart = 5.5
   tau_solve_accuracy_boost = 1.0
   timestep_boost = 1.0
   max_redshift = 50.0
 DarkEnergy: <DarkEnergyFluid>
   w = -1.0
   wa = 0.0
   cs2 = 1.0
   use_tabulated_w = False
 NonLinearModel: <Halofit>
   Min_kh_nonlinear = 0.005
   halofit_version = mead
   HMCode_A_baryon = 3.13
   HMCode_eta_baryon = 0.603
 Accuracy: <AccuracyParams>
   AccuracyBoost = 1.0
   lSampleBoost = 1.0
   lAccuracyBoost = 1.0
   AccuratePolarization = True
   AccurateBB = False
   AccurateReionization = True
   TimeStepBoost = 1.0
   BackgroundTimeStepBoost = 1.0
   IntTolBoost = 1.0
   SourcekAccuracyBoost = 1.0
   IntkAccuracyBoost = 1.0
   TransferkBoost = 1.0
   NonFlatIntAccuracyBoost = 1.0
   BessIntBoost = 1.0
   LensingBoost = 1.0
   NonlinSourceBoost = 1.0
   BesselBoost = 1.0
   LimberBoost = 1.0
   SourceLimberBoost = 1.0
   KmaxBoost = 1.0
   neutrino_q_boost = 1.0
 SourceTerms: <SourceTermParams>
   limber_windows = True
   limber_phi_lmin = 100
   counts_density = True
   counts_redshift = True
   counts_lensing = False
   counts_velocity = True
   counts_radial = False
   counts_timedelay = True
   counts_ISW = True
   counts_potential = True
   counts_evolve = False
   line_phot_dipole = False
   line_phot_quadrupole = False
   line_basic = True
   line_distortions = True
   line_extra = False
   line_reionization = False
   use_21cm_mK = True
 z_outputs = []
 scalar_initial_condition = initial_adiabatic
 InitialConditionVector = []
 OutputNormalization = 1
 Alens = 1.0
 MassiveNuMethod = Nu_best
 DoLateRadTruncation = True
 Evolve_baryon_cs = False
 Evolve_delta_xe = False
 Evolve_delta_Ts = False
 Do21cm = False
 transfer_21cm_cl = False
 Log_lvalues = False
 use_cl_spline_template = True
 SourceWindows = []
 CustomSources: <CustomSources>
   num_custom_sources = 0
   c_source_func = None
   custom_source_ell_scales = []
 Q0 = 0.0
 Qinf = 0.0
 DR0 = 0.0
 DRinf = 0.0
 s = 0.0
 k_c = 0.01
 Q1 = 1.5
 Q2 = 0.5
 Q3 = 1.5
 Q4 = 0.5
 D1 = 1.0
 D2 = 1.0
 D3 = 1.0
 D4 = 1.0
 z_div = 1.0
 z_TGR = 2.0
 z_tw = 0.05
 k_tw = 0.001
 t_dep = False
 R_func = False
 true_bin = True
 ISiTGR_BIN = True
 E11 = 0.0
 E22 = 0.0
 mu0 = 0.0
 Sigma0 = 0.0
 c1 = 1.0
 c2 = 1.0
 Lambda = 0.0
 mu1 = 0.0
 mu2 = 0.0
 mu3 = 0.0
 mu4 = 0.0
 eta1 = 0.0
 eta2 = 0.0
 eta3 = 0.0
 eta4 = 0.0
 Sigma1 = 0.0
 Sigma2 = 0.0
 Sigma3 = 0.0
 Sigma4 = 0.0
 ISiTGR_mueta = False
 ISiTGR_muSigma = False
 ISiTGR_QDR = False
 ISiTGR_BIN_mueta = False
 ISiTGR_BIN_muSigma = False
 GR = 1
In [23]:
# Getting CAMB results
results_GR = isitgr.get_results(pars_GR)
results_hybrid = isitgr.get_results(pars_binning_hybrid)
results_trad = isitgr.get_results(pars_binning_trad)
In [26]:
# Non-Linear spectra (Halofit)
#pars_GR.NonLinear = model.NonLinear_both
#results_GR.calc_power_spectra(pars_GR)
#kh_nonlin, z_nonlin, pk_nonlin = results_GR.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)

#Linear spectra
pars_GR.NonLinear = model.NonLinear_none
kh_GR, z_GR, pk_GR = results_GR.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)

#Linear spectra for MG (hybrid binning)
pars_binning_hybrid.NonLinear = model.NonLinear_none
kh_hybrid, z_hybrid, pk_hybrid = results_hybrid.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)

#Linear spectra for MG (traditional binning)
pars_binning_trad.NonLinear = model.NonLinear_none
kh_trad, z_trad, pk_trad = results_trad.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)
In [27]:
print(results_GR.get_sigma8())
print(results_hybrid.get_sigma8())
print(results_trad.get_sigma8())

[0.49106338 0.79068885]
[0.51066247 0.91396426]
[0.51089123 0.91549571]
In [28]:
ISiTGR_pk_hybrid = pk_hybrid[0,:]#*0.7 #normalizing P(k)
ISiTGR_pk_trad = pk_trad[0,:]#*0.7 #normalizing P(k)
Cris_PK = pk_GR[0,:] * 1.35

plt.figure(figsize=(10,10))
plt.loglog(kh_GR, Cris_PK, color='k', label='$\Lambda CDM$')
plt.loglog(kh_hybrid, ISiTGR_pk_hybrid, color='r', linestyle = 'dashdot', label='Hybrid Binning')
plt.loglog(kh_trad, ISiTGR_pk_trad, color='b', linestyle = 'dashed',label='Traditional Binning')

plt.xlim([0.001,1])
plt.legend(loc='lower left', fontsize=16)
plt.ylabel('$P_\delta (k)$', fontsize=16)
plt.xlabel('$k$', fontsize=16)
plt.title('obtained using ISiTGR-python wrapper')
Out[28]:
Text(0.5, 1.0, 'obtained using ISiTGR-python wrapper')
In [ ]:
In [ ]: