Initializing session

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from classy_sz import Class
from scipy.optimize import fsolve
from scipy.interpolate import interp1d
import math

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

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


cosmology = {
                   # LambdaCDM parameters
                   'H0':67.556,
                   'omega_b':0.022032,
                   'omega_cdm':0.12038,
                   # 'sigma8': 0.831,
                   'ln10^{10}A_s': 3.047,
                   'n_s': 0.9665,
                   'tau_reio':0.0925,
}


import os
path_to_class_sz = os.environ['PATH_TO_CLASS_SZ_DATA']

Galaxy kernels

z , nz_lenses = np.loadtxt(path_to_class_sz+'/class_sz/class-sz/class_sz_auxiliary_files/includes/nz_lens_bin1.txt',unpack=True)
plt.plot(z,nz_lenses,label='lenses')
print(np.trapz(nz_lenses,z))
z , nz_sources = np.loadtxt(path_to_class_sz+'/class_sz/class-sz/class_sz_auxiliary_files/includes/nz_source_normalized_bin4.txt',unpack=True)
plt.plot(z,nz_sources,label='sources (normalized)',ls='--')

print(np.trapz(nz_sources,z))
plt.legend()

0.9999999999999999
1.0000000000000002

<matplotlib.legend.Legend at 0x34d179250>

Compute

%%time

M = Class()
M.set(cosmology)
M.set({
'output':'gallens_lens_1h,gallens_lens_2h',
'ell_max' : 5e3,
'ell_min' : 2,
'dlogell':0.1,


'ndim_redshifts':80,

'mass_function':'T08M200c',


'M_min':1e11, # Msun/h
'M_max':1e15, # Msun/h
# 'mass_epsrel': 1e-3,
# 'mass_epsabs': 1e-40,

'z_min':1e-5,
'z_max': 2.,
# 'redshift_epsrel': 0.5e-3,
# 'redshift_epsabs': 1e-40,

'delta_for_galaxies':'200c',
'delta_for_matter_density':'200c',

'concentration_parameter':'fixed',


'Delta_z_source':0.00,

'galaxy_sample' : 'custom',
'full_path_to_source_dndz_gal' : path_to_class_sz+'/class_sz/class-sz/class_sz_auxiliary_files/includes/nz_source_normalized_bin4.txt', # source galaxies

'N_samp_fftw':1024, #precision parameter for the bessel transform to theta space
'l_min_samp_fftw' : 1e-8,
'l_max_samp_fftw' : 1e8,

'hm_consistency' : 1,

})
M.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 25.5 s, sys: 2.17 s, total: 27.7 s
Wall time: 3.05 s

ell =  np.asarray(M.cl_kg_k()['ell'])
cl1h = M.cl_kg_k()['1h']
cl2h = M.cl_kg_k()['2h']

Plot

label_size = 17
title_size = 18
legend_size = 13
handle_length = 1.5
fig, (ax1) = plt.subplots(1,1,figsize=(6,4))

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)
ax.grid( visible=True, which="both", alpha=0.2, linestyle='--')
ax.set_xlabel(r'$\ell$',size=title_size)
ax.set_ylabel(r'$C_\ell^{\kappa\gamma}$',size=title_size)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim(1e-11,1.5e-7)
ax.set_xlim(2.,5e3)
fac = ell*(ell+1.)/2./np.pi
ax.plot(ell,cl1h/fac,label='CLASS_SZ 1h',c='b',ls='-',lw=0.7)
ax.plot(ell,cl2h/fac,label='CLASS_SZ 2h',c='b',ls='--',lw=0.7)
ax.plot(ell,cl1h/fac+cl2h/fac,label='CLASS_SZ 1+2h',c='r')


ax.legend(loc=1,ncol = 1,frameon=False,fontsize=11)

fig.tight_layout()