

In [1]:

    
import camb
pars = camb.CAMBparams()
from scipy.special import jn, jn_zeros
from camb import model, initialpower
from scipy.interpolate import interp1d
from hankel_transform import *
from astropy.cosmology import Planck13 #use Planck15 if you can
import astropy.units as u
#rc('text', usetex=False)









    



---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-0f9c0e9836b1> in <module>()
----> 1 import camb
      2 pars = camb.CAMBparams()
      3 from scipy.special import jn, jn_zeros
      4 from camb import model, initialpower
      5 from scipy.interpolate import interp1d

ImportError: No module named camb



In [ ]:

    
from matplotlib.pyplot import *
%matplotlib inline



In [ ]:

    
cosmo=Planck13.clone(H0=100)#we want h=1 for some comparisons



In [ ]:

    
#Set up a new set of parameters for CAMB
kmax=20#30
kmin=.8e-3
k_smooth=1
nk=5000
rmin=.1
rmax=110
non_linear=1

pars = camb.CAMBparams()
#This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency
pars.set_cosmology(H0=67.27, ombh2=0.022250, omch2=0.119800, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(ns=0.965, r=0,As =2.14e-09)
pars.set_for_lmax(2500, lens_potential_accuracy=0)
zb=[.27]#[0.0]
pars.set_matter_power(redshifts=zb,kmax=kmax);
if non_linear==1:
    pars.NonLinear = model.NonLinear_both
else:
    pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk =results.get_matter_power_spectrum(minkh=kmin, maxkh=kmax, npoints =nk)



In [ ]:

    
#Setting up the Hankel Transform
#This part is slower. But only needs to be run once. 
#If you only need wgg, set j_nu=[0]. For wg+ (or \Delta\Sigma) use j_nu=[2]
%time HT=hankel_transform(rmin=rmin,rmax=rmax,kmax=kmax,j_nu=[0],n_zeros=28000,kmin=kmin)
#HT=hankel_transform(rmin=1,rmax=rmax,kmax=1,j_nu=[0,2],n_zeros=2800,kmin=1.e-2)#quick test... inaccurate



In [ ]:

    
pk_taper=HT.taper(k=kh,pk=pk[0],large_k_lower=5,large_k_upper=kmax,low_k_lower=kmin,
                     low_k_upper=kmin*2)
#need atleast k=10 and k=1.e-3 to get decent wgg for 1-100 Mpc/h range.



In [ ]:

    
b_g=1



In [ ]:

    
r_gg,wgg=HT.projected_correlation(k_pk=kh,pk=pk[0]*b_g**2,j_nu=0)
%time r_gg,wgg_taper=HT.projected_correlation(k_pk=kh,pk=pk_taper*b_g**2,j_nu=0)

#Internal test case wgg_test_fact=1./np.pi**2 wgg_test=np.genfromtxt('../data/planck15_wgg_rsd0_kz0_s0_k10_z0.27.dat',names=('rp','wgg')) wgp_test=np.genfromtxt('../data/planck15_wgp_rsd0_kz0_s0_k10_z0.27.dat',names=('rp','wgp'))

plot(kh,pk[0]) plot(kh,pk_taper) xscale('log') yscale('log') xlim(1.e-4,30) ylim(1,5.e4)



In [ ]:

    
plot(r_gg,r_gg*wgg,label='No tapering')
plot(r_gg,r_gg*wgg_taper,label='tapering')
#plot(wgg_test['rp'],wgg_test['rp']*wgg_test['wgg']*wgg_test_fact)
xscale('log')
#yscale('log')
#ylim(1.e-2,1.e3)
xlim(1,110)
xlabel('$r_p$ [Mpc/h]')
ylabel('$w_{gg}$ [Mpc/h]')
ylim(0,72)
legend()



In [ ]:

    
plot(r_gg,wgg-wgg_taper,'r')
xscale('log')
xlabel('$r_p$ [Mpc/h]')

w_intp=interp1d(r,w,bounds_error=False,fill_value=np.nan) plot(wgg_test['rp'],w_intp(wgg_test['rp'])/(wgg_test['wgg']*wgg_test_fact),'r--') w_intp=interp1d(r,w2,bounds_error=False,fill_value=np.nan) plot(wgg_test['rp'],w_intp(wgg_test['rp'])/(wgg_test['wgg']*wgg_test_fact)) xscale('log') ylim(.98,1.02)

g-g Lensing covariance



In [ ]:

    
%time HT=hankel_transform(rmin=rmin,rmax=rmax,kmax=kmax,j_nu=[2],n_zeros=60000,kmin=kmin)



In [ ]:

    
sigma_crit=4.7e3



In [ ]:

    
sigma_e=.36
sigma_gamma=sigma_e/1.7
n_s=8
shape_noise=sigma_crit**2*sigma_gamma**2/n_s#*cosmo.H_z(z=0.27)/cosmo.c



In [ ]:

    
rho=0.04180594698596614 #units of Msun/pc^2/mpc, using h=1



In [ ]:

    
Dchi2=900 #integral over product of window functions
p_kappa_kappa=pk[0]*(rho)**2*Dchi2



In [ ]:

    
n_g=5e-4#3.e-4
g_shot_noise=1./n_g
b_g=1.838
p_g=b_g**2*pk[0]
#p_gk=b_g*pk[0]*rho



In [ ]:

    
r_bins=np.logspace(-1,1.6,19)



In [ ]:

    
area=10000
area_comoving=area*(np.pi/180)**2*cosmo.comoving_distance(z=.27)**2
L_W=500
vol=area_comoving*L_W
vol=vol.value

from pearce.mocks.kittens import DarkSky cat = DarkSky(264, system = 'sherlock')



In [ ]:

    
pars = camb.CAMBparams()
#This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency

# replaced with DarkSky values
pars.set_cosmology(H0=70.4, ombh2=0.0231566, omch2=0.12293, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(ns=0.967, r=0,As =2.18e-09)
pars.set_for_lmax(2500, lens_potential_accuracy=0)
zb=[0.0]#[0.0]#
pars.set_matter_power(redshifts=zb,kmax=kmax);
if non_linear==1:
    pars.NonLinear = model.NonLinear_both
else:
    pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk =results.get_matter_power_spectrum(minkh=kmin, maxkh=kmax, npoints =nk)

#real_p_g = np.load('/home/users/swmclau2/Git/pearce/bin/covmat/galaxy_pk.npy') real_p_g = np.load('/home/users/swmclau2/Git/pearce/notebooks/p_g.npy')

real_p_g = real_p_g[1:]



In [ ]:

    
taper_kw=dict({'large_k_lower':10,'large_k_upper':kmax,'low_k_lower':kmin,'low_k_upper':kmin*1.2})



In [ ]:

    
r,cov_ggkk=HT.projected_covariance(k_pk=kh,pk1=p_g+g_shot_noise,pk2=np.zeros_like(p_kappa_kappa)+shape_noise,j_nu=2,taper=True,**taper_kw)
r_re,cov_ggkk_re=HT.bin_cov(r=r,cov=cov_ggkk,r_bins=r_bins)
corr=HT.corr_matrix(cov=cov_ggkk_re)

cov_ggkk_re/=vol # Don't forget!

plot(kh, pk.squeeze()) plot(real_p_g[:,0], real_p_g[:,1]) loglog();

np.sqrt(real_p_g[np.argmin(np.abs(real_p_g[:,0] - 0.1) ),1]/pk.squeeze()[np.argmin(np.abs(kh - 0.1) )])

r,cov_ggkk=HT.projected_covariance(k_pk=real_p_g[:,0],pk1=real_p_g[:,1]+g_shot_noise,pk2=np.zeros_like(real_p_g[:,0]) + shape_noise,j_nu=2,taper=True,**taper_kw) r_re,cov_ggkk_re=HT.bin_cov(r=r,cov=cov_ggkk,r_bins=r_bins) corr=HT.corr_matrix(cov=cov_ggkk_re)



In [ ]:

    
plot(r_re, np.sqrt(np.diag(cov_ggkk_re)) )
loglog();



In [ ]:

    
print np.sqrt(np.diag(cov_ggkk_re))



In [ ]:

    
imshow(corr)



In [ ]:

    
np.save('shape_noise.npy', cov_ggkk_re)



In [76]:

    
np.min(r_re), np.min(r)









    Out[76]:





(0.11809247432396976, 0.1)



In [38]:

    
from matplotlib.colors import LogNorm



In [39]:

    
pcolor(r_re,r_re,corr,norm=LogNorm())
colorbar()
xscale('log')
yscale('log')
xlim(HT.rmin,HT.rmax)
ylim(HT.rmin,HT.rmax)









    Out[39]:





(0.1, 110)



In [40]:

    
r,cov_gkgk=HT.projected_covariance(k_pk=kh,pk1=p_gk,pk2=p_gk,
                               kmax=100,rmin=.8,j_nu=2,rmax=110,n_zeros=3500)#return_Jrr=True,Jrr=Jrr
cov_gkgk*=Dchi2
r_re,cov_gkgk_re=HT.bin_cov(r=r,cov=cov_gkgk,r_bins=r_bins)
corr=HT.corr_matrix(cov=cov_gkgk_re)



In [41]:

    
cov_final=(cov_ggkk_re+cov_gkgk_re)/vol
corr=HT.corr_matrix(cov=cov_final)
errors=HT.diagonal_err(cov=cov_final)



In [42]:

    
pcolor(r_re,r_re,corr,norm=LogNorm())
colorbar()
xscale('log')
yscale('log')
xlim(HT.rmin,HT.rmax)
ylim(HT.rmin,HT.rmax)









    Out[42]:





(0.1, 110)

error_qpm=np.genfromtxt('../data/qpm_DS_error_th.dat',names=True)



In [43]:

    
plot(r_re,r_re*errors,label='new')
#plot(error_qpm['rp'],error_qpm['rp']*error_qpm['DS_err_gR'])
xscale('log')
yscale('log')
legend()
xlim(1,100)









    Out[43]:





(1, 100)



In [ ]:



In [ ]: