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