In [1]:

    

from power_spectra import * #look into this file. You need camb or class to get power spectra (there are wrapper functions for both). 
                            #The file imports both class and camb. Comment out relevant lines (at the top) if you don't have one of them.
from scipy.special import jn, jn_zeros
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)
import scipy
from scipy.interpolate import interp1d


    

In [2]:

    

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


    

In [3]:

    

#set parameters for power spectra computations.
#play around with kmax,kmin to find optimal values for your case (based on rmin,rmax). 
#Remember very high kmax or very low kmin will require finer grid, slowing down the initial setup of the hankel transform. 
#Subsequent calculations of correlation functions are fast, so speed should not be major issue when running fits.
#need atleast k=10 and k=1.e-3 to get decent wgg for 1-100 Mpc/h range.
kmax=30
kmin=.8e-3
k_smooth=1
nk=5000
rmin=1
rmax=100
non_linear=1

r_bins=np.logspace(-1,np.log10(100),21) #np.logspace(0,2,11)


    

In [4]:

    

cosmo_fid=dict({'h':cosmo.h,'Omb':cosmo.Ob0,'Omd':cosmo.Om0-cosmo.Ob0,'s8':0.817,'Om':cosmo.Om0,
                'As':2.12e-09,'mnu':cosmo.m_nu[-1].value,'Omk':cosmo.Ok0,'tau':0.06,'ns':0.965,'w':-1,'wa':0})
pk_params={'non_linear':non_linear,'kmax':kmax,'kmin':kmin,'nk':nk}


    

In [5]:

    

PS=Power_Spectra(cosmo_params=cosmo_fid,pk_params=pk_params)


    

In [6]:

    

sc1=PS.sigma_crit(zl=np.atleast_1d(0.27),zs=np.atleast_1d(0.45),cosmo_h=cosmo_h)
Oms=np.linspace(cosmo.Om0*.7,cosmo.Om0*1.3,11)
scs=np.zeros_like(Oms)
i=0
for om in Oms:
    cosmo_i=cosmo_h.clone(Om0=om)
    scs[i]=PS.sigma_crit(zl=np.atleast_1d(0.27),zs=np.atleast_1d(0.45),cosmo_h=cosmo_i).value[0,0]
    i+=1


    

In [7]:

    

plot(Oms/Oms[5],scs/scs[5])
plot(Oms/Oms[5],(Oms/Oms[5])*.12+.88,'--')
# loglog()


    

    
Out[7]:





[<matplotlib.lines.Line2D at 0x7fd958b28400>]






    




/Users/Deep/anaconda/envs/py36/lib/python3.6/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))






    

In [8]:

    

z_min=0.16
z_max=0.36
z_mean=0.27


    

In [9]:

    

# pk,kh =PS.class_pk(z=[z_mean])
pk,kh =PS.camb_pk(z=[z_mean])


    

In [10]:

    

#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,2,4],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


    

    




j-nu= 0  not enough zeros to cover kmin, increasing by  1000  to 29000
nr: 948
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 29000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 30000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 31000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 32000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 33000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 34000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 35000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 36000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 37000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 38000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 39000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 40000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 41000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 42000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 43000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 44000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 45000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 46000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 47000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 48000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 49000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 50000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 51000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 52000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 53000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 54000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 55000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 56000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 57000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 58000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 59000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 60000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 61000
j-nu= 2  not enough zeros to cover kmin, increasing by  1000  to 62000
nr: 948
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 29000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 30000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 31000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 32000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 33000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 34000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 35000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 36000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 37000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 38000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 39000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 40000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 41000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 42000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 43000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 44000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 45000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 46000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 47000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 48000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 49000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 50000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 51000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 52000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 53000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 54000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 55000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 56000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 57000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 58000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 59000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 60000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 61000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 62000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 63000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 64000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 65000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 66000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 67000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 68000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 69000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 70000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 71000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 72000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 73000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 74000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 75000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 76000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 77000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 78000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 79000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 80000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 81000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 82000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 83000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 84000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 85000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 86000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 87000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 88000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 89000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 90000
j-nu= 4  not enough zeros to cover kmin, increasing by  1000  to 91000
nr: 948
CPU times: user 53.3 s, sys: 906 ms, total: 54.2 s
Wall time: 54.3 s

In [11]:

    

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)
#tapering helps to reduce ringing in the final correlation functions


    

In [12]:

    

from scipy.integrate import quad as scipy_int1d


    

In [13]:

    

cosmo.H(z=[0]).value


    

    
Out[13]:





array([67.77])

In [14]:

    

def DZ_int(z=[0],cosmo=None,rtol=1.e-4,tol=1.e-5): #linear growth factor.. full integral.. eq 63 in Lahav and suto
        Ez_func=cosmo.efunc
        def intf(z):
            return (1.+z)/(cosmo.H(z=z).value**3)
        dz=np.zeros_like(z,dtype='float32')
        inf=np.inf
        j=0
        for i in z:
            dz[j]+=cosmo.H(i).value*scipy_int1d(intf,i,inf,epsrel=rtol,epsabs=tol)[0]
            j=j+1
        dz=dz*2.5*cosmo.Om0*cosmo.H0**2
        return dz/dz[0] #check for normalization


    

In [15]:

    

b_g=1.8
A_I=4.5
C1_rhoC=0.0134
Om=cosmo.Om0
Dz=DZ_int(z=np.append([0],z_mean),cosmo=cosmo)# Dz is normalized to 1 at lowest z. 
                                    #Hence appending 0 in front to get Dz normalized to 1 at z=0
Dz=Dz[1]/Dz[0]

#prefactors in front of wmm or wm+.
#if doing fits, set b_g and AI to be 1 here and then in the fitting function 
# you can simply multiply the correlation function with the correct values.
wgg_f=b_g**2
wgp_f=b_g*A_I*C1_rhoC*Om/Dz
wpp_f=(A_I*C1_rhoC*Om/Dz)**2


    

In [16]:

    

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

r_gp,wgp=HT.projected_correlation(k_pk=kh,pk=pk[0]*wgp_f,j_nu=2)

r_pp,wpp4=HT.projected_correlation(k_pk=kh,pk=pk[0]*wpp_f,j_nu=4)
r_pp0,wpp0=HT.projected_correlation(k_pk=kh,pk=pk[0]*wpp_f,j_nu=0)

wpp0_intp=interp1d(r_pp0,wpp0,bounds_error=False,fill_value=np.nan)
wpp=wpp4+wpp0_intp(r_pp)


    

    




CPU times: user 69 ms, sys: 1.17 ms, total: 70.2 ms
Wall time: 8.9 ms

In [17]:

    

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]')

legend()


    

    
Out[17]:





<matplotlib.legend.Legend at 0x7fd9486738d0>






    

In [18]:

    

plot(r_gp,wgp,label='No tapering')
xscale('log')
yscale('log')
xlim(1,110)
xlabel('$r_p$ [Mpc/h]')
ylabel('$w_{g+}$ [Mpc/h]')
ylim(1.e-2,20)
legend()


    

    
Out[18]:





<matplotlib.legend.Legend at 0x7fd978bb0048>






    

In [19]:

    

plot(r_pp,r_pp*wpp,label='No tapering')
xscale('log')
xlim(1,110)
xlabel('$r_p$ [Mpc/h]')
ylabel('$w_{++}$ [Mpc/h]')
# ylim(0,72)
legend()


    

    
Out[19]:





<matplotlib.legend.Legend at 0x7fd958bd0940>






    

In [20]:

    

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


    

    
Out[20]:





Text(0.5, 0, '$r_p$ [Mpc/h]')






    

In [21]:

    

n_g=3.e-4 #number density.. in units of Mpc ^ -3, same as pk
g_shot_noise=1./n_g
b_g=1.8
p_g=pk[0]*wgg_f


    

In [22]:

    

sigma_e=.36
sigma_gamma=sigma_e/1.7
n_s=n_g #shape galaxy number density 
shape_noise=sigma_gamma**2/n_s#*cosmo.H_z(z=0.27)/cosmo.c
p_gk=pk[0]*wgp_f
p_kk=pk[0]*wpp_f


    

In [23]:

    

#lensing contribution. For sims, you can set p_kappa to 0
l=np.unique(np.int32(kh*cosmo.comoving_transverse_distance(z_mean).value))
l,cl=PS.kappa_cl(zs1=[z_mean],zs2=[z_mean],p_zs1=[1],p_zs2=[1],pk_func=PS.camb_pk,zl_max=z_mean,l=l) #set to class_pk if using class
p_kappa=cl*cosmo.comoving_transverse_distance(z_mean)**3 / 2 #this is an approximation
k_kappa=l/cosmo.comoving_transverse_distance(z_mean)
p_kappa_int=interp1d(k_kappa,p_kappa,bounds_error=False, fill_value=0)
p_kappa=p_kappa_int(kh)


    

    




Note: redshifts have been re-sorted (earliest first)






    




/Users/Deep/BOX/Box Sync/repos/Public-code/Hankel_transform/power_spectra.py:242: RuntimeWarning: divide by zero encountered in true_divide
  pk_int=interp1d(lz,pk[i]/DC_i**2,bounds_error=False,fill_value=0)
/Users/Deep/anaconda/envs/py36/lib/python3.6/site-packages/astropy/units/quantity.py:463: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/Deep/BOX/Box Sync/repos/Public-code/Hankel_transform/power_spectra.py:274: RuntimeWarning: divide by zero encountered in true_divide
  f=(l+0.5)**2/(l*(l+1.)) #correction from Kilbinger+ 2017

In [24]:

    

plot(kh,p_kappa,'--')
plot(kh,p_kk)
loglog()


    

    
Out[24]:





[]






    

In [25]:

    

area=8000 #sq. degrees
area_comoving=area*(np.pi/180)**2*cosmo.comoving_distance(z=z_mean)**2
L_W=cosmo.comoving_distance(z=z_max)-cosmo.comoving_distance(z=z_min)
vol=area_comoving*L_W
vol=vol.value
area_comoving=area_comoving.value

#for sims, use the box dimensions to get area and volume.


    

In [26]:

    

p_g_cov=p_g+g_shot_noise
p_kk_cov=p_kk+p_kappa+shape_noise
p_gk_cov=p_gk


    

In [27]:

    

r,cov_ggkk=HT.projected_covariance(k_pk=kh,pk1=p_g_cov,pk2=p_kk_cov,j_nu=2)
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)


    

    




/Users/Deep/BOX/Box Sync/repos/Public-code/Hankel_transform/hankel_transform.py:130: RuntimeWarning: invalid value encountered in double_scalars
  corr[i][j]=cov[i][j]/np.sqrt(cov[i][i]*cov[j][j])

In [28]:

    

r,cov_gkgk=HT.projected_covariance(k_pk=kh,pk1=p_gk_cov,pk2=p_gk_cov,j_nu=2)
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 [29]:

    

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


    

In [30]:

    

lowz_dat=np.genfromtxt('./test_Dat/lowz_full.dat',names=True)
lowz_cov=np.genfromtxt('./test_Dat//lowz_wgp_cov.dat')
lowz_corr=HT.corr_matrix(cov=lowz_cov)


    

In [31]:

    

pcolor(r_re,r_re,corr,vmin=-1,vmax=1,cmap='seismic')
colorbar()
xscale('log')
yscale('log')
xlim(HT.rmin,HT.rmax)
ylim(HT.rmin,HT.rmax)
# colorbar()
show()
pcolor(lowz_dat['rp'],lowz_dat['rp'],corr,vmin=-1,vmax=1,cmap='seismic')
colorbar()
xscale('log')
yscale('log')
xlim(HT.rmin,HT.rmax)
ylim(HT.rmin,HT.rmax)
# colorbar()


    

    













    
Out[31]:





(1, 100)






    

In [32]:

    

plot(r_gp,wgp)
plot(r_re,errors,'--')
xscale('log')
yscale('log')
xlim(1,110)


    

    
Out[32]:





(1, 110)






    

In [34]:

    

plot(r_gp,wgp)
plot(lowz_dat['rp'],lowz_dat['wgp'])
plot(lowz_dat['rp'],np.sqrt(np.diag(lowz_cov)))
plot(r_re,errors)
xscale('log')
yscale('log')
xlim(1,110)
ylim(9.e-3,2)


    

    
Out[34]:





(0.009, 2)






    

In [ ]:

    




    

In [ ]: