In [1]:
from hankel_transform import *
from power_spectra import *
from astropy.cosmology import Planck15 as cosmo
In [2]:
PS=Power_Spectra()
In [3]:
lmin=2
lmax=20000
d2r=np.pi/180.#degree to rad
a2r=1./60*d2r #arcmin_to_rad
theta_min=1.*a2r
theta_max=600.*a2r
theta_bins=np.logspace(0,np.log10(600),21)#bin_edges
In [4]:
HT=hankel_transform(rmin=theta_min,rmax=theta_max,kmax=lmax,kmin=lmin,j_nu=[0],
n_zeros=16000,prune_r=2)
lmax=int(max(HT.k[0])+1)
In [5]:
l=np.arange(lmax,dtype='int')
In [6]:
#cross power spectra. Cl_XY in notes
l,kappa_cl=PS.kappa_cl(n_zl=20,log_zl=False,zl_min=1.e-2,zl_max=1,zs1=[1],p_zs1=[1],zs2=[1.5],p_zs2=[1],l=l)
In [7]:
#auto terms.. Cl_XX and Cl_YY in notes
l,kappa1_cl=PS.kappa_cl(n_zl=20,log_zl=False,zl_min=1.e-2,zl_max=1,zs1=[1],p_zs1=[1],zs2=[1],p_zs2=[1],l=l)
l,kappa2_cl=PS.kappa_cl(n_zl=20,log_zl=False,zl_min=1.e-2,zl_max=1,zs1=[1.5],p_zs1=[1],zs2=[1.5],p_zs2=[1],l=l)
In [8]:
plot(l,kappa_cl)
xscale('log')
yscale('log')
In [9]:
taper_kw=dict({'large_k_lower':lmax/2.,'large_k_upper':lmax,'low_k_lower':lmin,
'low_k_upper':lmin*1.2})
In [10]:
n_s=1 #per sq arcmin
n_s=n_s/a2r**2
sigma_e=0.36
shape_noise=sigma_e**2/n_s
N_cl1=np.ones_like(kappa_cl)*shape_noise
#N_cl1*=100
N_cl2=N_cl1
In [11]:
shape_noise
Out[11]:
In [12]:
theta,kappa_corr=HT.projected_correlation(k_pk=l,pk=kappa_cl,j_nu=0,taper=True,**taper_kw)
In [13]:
plot(theta/a2r,kappa_corr)
xscale('log')
yscale('log')
In [14]:
#Cl_XY**2 term in covariance
theta,cov_cross=HT.projected_covariance(k_pk=l,pk1=kappa_cl,pk2=kappa_cl,j_nu=0,taper=True,**taper_kw)
In [ ]:
#Noise covariance.. Not needed, just to check code. It should go as 1/theta if everything is working. There
#might be some resolution effects if we are not using very fine grid in HT.
theta,cov_cross=HT.projected_covariance(k_pk=l,pk1=kappa_cl,pk2=kappa_cl,j_nu=0,taper=True,**taper_kw)
In [15]:
#binning the covariance.
theta_b,cov_cross_b=HT.bin_cov(r=theta/a2r,cov=cov_cross,r_bins=theta_bins)
theta_b,N_cov_b=HT.bin_cov(r=theta/a2r,cov=N_cov,r_bins=theta_bins)
In [16]:
plot(theta_b,np.sqrt(np.diagonal(cov_cross_b)))
plot(theta_b,np.sqrt(np.diagonal(N_cov_b)),label='N')
xscale('log')
yscale('log')
xlabel(r'$\theta$ (arc min)')
ylabel(r'$\delta\kappa$')
legend()
Out[16]:
In [17]:
#Cl_XX*Cl_YY term in covariance
theta,cov_auto=HT.projected_covariance(k_pk=l,pk1=kappa1_cl+N_cl1,pk2=kappa2_cl+N_cl2,j_nu=0,
taper=True,**taper_kw)
In [18]:
cov_full=cov_auto+cov_cross
theta_b,cov_full_b=HT.bin_cov(r=theta/a2r,cov=cov_full,r_bins=theta_bins)
In [19]:
kappa_int=interp1d(theta/a2r,kappa_corr)
kappa_b=kappa_int(theta_b)
In [21]:
np.diagonal(cov_full_b)
Out[21]:
In [24]:
errorbar(theta_b,kappa_b,np.sqrt(np.diagonal(cov_full_b)),fmt='bo')
xscale('log')
yscale('log')
xlabel(r'$\theta$ (arc min)')
ylabel(r'$\kappa\kappa$')
legend()
In [26]:
plot(theta_b,np.sqrt(np.diagonal(cov_full_b)),'b-')
xscale('log')
yscale('log')
xlabel(r'$\theta$ (arc min)')
ylabel(r'$\delta \kappa\kappa$')
legend()
In [27]:
#Cl_XX Cl_YY Cl_XY term in notes
theta,skew_auto_cross=HT.skewness(k_pk=l,pk1=kappa1_cl+N_cl1,pk2=kappa2_cl+N_cl2,pk3=kappa_cl,
j_nu=0,taper=True,**taper_kw)
#this is annoyingly slow right now, because of einsum statement. This isn't really einsum's fault though.
#Calculation is N^3 here vs N^2 in covariance.
In [28]:
#Cl_XY**3 term in notes
theta,skew_cross=HT.skewness(k_pk=l,pk1=kappa_cl,pk2=kappa_cl,pk3=kappa_cl,j_nu=0,taper=True,**taper_kw)
In [29]:
skew_full=(skew_auto_cross)*6+(skew_cross)*2
In [30]:
#bin 3rd moment
theta_b,skew_full_b=HT.bin_mat(r=theta/a2r,mat=skew_full,r_bins=theta_bins)
In [31]:
def skew_corr(skew=[]):
ndim=len(skew)
skew_corr=np.zeros_like(skew)
for i in np.arange(ndim):
for j in np.arange(ndim):
for k in np.arange(ndim):
skew_corr[i][j][k]=skew[i][j][k]/(skew[i][i][i]*skew[j][j][j]*skew[k][k][k])**(1./3.)
return skew_corr
In [32]:
#correlation matrix like term, but for 3rd moment.
scorr=skew_corr(skew_auto_cross)
scorr_b=skew_corr(skew_full_b)
In [33]:
# pretty plot, looking at scorr along one axis.
pcolor(theta/a2r,theta/a2r,(scorr[0]),vmax=1,vmin=-1,cmap='seismic')
xscale('log')
yscale('log')
In [34]:
plot(theta_b,np.diagonal(np.diagonal(skew_full_b)))
plot(theta/a2r,np.diagonal(np.diagonal(skew_full)))
plot(theta_b,np.diagonal(cov_full_b**1.5))
xscale('log')
yscale('log')
In [35]:
a=plot(theta/a2r,np.diagonal(np.diagonal(skew_cross))/np.diagonal(cov_cross)**1.5)
a=plot(theta/a2r,np.diagonal(np.diagonal(skew_auto_cross))/np.diagonal(cov_cross)**1.5)
xscale('log')
In [36]:
a=plot(theta/a2r,np.diagonal(np.diagonal(skew_full))/np.diagonal(cov_full)**1.5)
a=plot(theta_b,np.diagonal(np.diagonal(skew_full_b))/np.diagonal(cov_full_b)**1.5)
xscale('log')
In [ ]: