Correlation function of DR72 SDSS VAGC Catalog

First import all the modules such as healpy and astropy needed for analyzing the structure



In [1]:

    
import healpix_util as hu
import astropy as ap
import numpy as np
from astropy.io import fits
from astropy.table import Table
import astropy.io.ascii as ascii
from astropy.io import fits
from astropy.constants import c
import matplotlib.pyplot as plt
import math as m
from math import pi
#from scipy.constants import c
import scipy.special as sp
from astroML.decorators import pickle_results
from scipy import integrate
import warnings
from sklearn.neighbors import BallTree
import pickle
import multiprocessing as mp
import time
from aptestmetricdt import *
from aptestmetricdz import *
from scipy.spatial import distance as d
from apcat import *
from progressbar import *
from tqdm import *
from functools import partial
import pymangle
from scipy.optimize import curve_fit
#from astroML.datasets import fetch_sdss_specgals
#from astroML.correlation import bootstrap_two_point_angular
%matplotlib inline

Co-moving distance between any two sources at redshift $z$ separated by $\Delta z$ and $\Delta\theta$ in the redshift space is

$s(\Delta z, z\Delta\theta)\approxeq \frac{c}{H_0}x(z)\sqrt{(\Delta z)^2+y^2(z)(z\Delta\theta)^2}$

Here $x(z)$ and $y(z)$ are model independent parameters. $y(z)$ characterizes Alcock-Paczynski ratio.

$x(z)$ and $y(z)$ functions for $\Lambda CDM$ and Linear coasting models are

For $\Lambda CDM$

$x(z)=\frac{1}{\sqrt{\Omega_m(1+z)^3+\Omega_m}}$

$y(z)=\frac{1}{z}\int_0^z dx\sqrt{\frac{\Omega_m(1+z)^3+\Omega_\Lambda}{\Omega_m(1+x)^3+\Omega_\Lambda}}$

For Linear coasting in a flat universe

$x(z)=\frac{1}{E(z)}=\frac{1}{1+z}$

$y(z)=\frac{1+z}{z}\ln(1+z)$

For Milne Universe (open empty Universe that is also linearly coasting)

$x(z)=\frac{1}{E(z)}=\frac{1}{1+z}$

$y(z)=1+\frac{z}{2}$

$s_{\perp}=\frac{c}{H_0}x(z)y(z)z\Delta \theta$ $s_{\parallel}=\frac{c}{H_0}x(z)\Delta z$

Read the data file (taken from http://cosmo.nyu.edu/~eak306/SDSS-LRG.html ) converted to ascii with comoving distance etc. in V01 reading from pkl files for faster read



In [2]:

    
with open('dd2ddr72v06cdist.pkl') as f:
    dd2d=pickle.load(f)



In [3]:

    
with open('rr2ddr72v06cdist.pkl') as f:
    rr2d=pickle.load(f)



In [4]:

    
with open('dr2ddr72v06cdist.pkl') as f:
    dr2d=pickle.load(f)



In [5]:

    
corrells=(dd2d+rr2d-2.0*dr2d)/rr2d



In [6]:

    
corrells









    Out[6]:





array([[ 2.08273679, -0.06591932, -1.14525488, -1.49923222, -1.74361248,
        -1.81418446, -1.90115591, -1.94329043, -1.93130228, -1.96064259],
       [ 0.20180514, -0.94522815, -1.34635828, -1.58844515, -1.7528935 ,
        -1.84438752, -1.86257547, -1.9439332 , -1.99009407, -1.94829146],
       [-1.09749238, -1.48234452, -1.64234357, -1.75344742, -1.83019763,
        -1.88547473, -1.92958533, -1.940775  , -1.95535063, -1.9629588 ],
       [-1.62892483, -1.70176401, -1.76263555, -1.84966414, -1.87139241,
        -1.91879916, -1.90796703, -1.96632612, -1.95794984, -1.94585536],
       [-1.71551518, -1.91806287, -1.86999609, -1.90800839, -1.92515012,
        -1.93386969, -1.96487487, -1.95405297, -1.97032435, -1.96087146],
       [-1.81904322, -1.94790171, -2.        , -1.96932249, -1.91842475,
        -1.98101775, -1.98521938, -1.97618539, -1.95892852, -2.01410874],
       [-1.8138236 , -1.90877811, -1.93966404, -1.98406259, -1.95253807,
        -1.95867436, -1.96679456, -2.01126218, -1.99885932, -1.99369292],
       [-1.9927308 , -1.97658099, -1.91586737, -1.96813381, -1.98416011,
        -1.98639893, -1.99997964, -1.99746813, -1.98126505, -1.97767361],
       [-1.96406473, -1.92609437, -1.99171257, -1.94687523, -1.9823819 ,
        -1.97786533, -2.02323228, -1.98559114, -1.99172681, -2.00881064],
       [-1.92470024, -1.99537152, -1.99419303, -1.97261023, -1.9673841 ,
        -1.98801505, -2.02412735, -2.01394008, -2.01956817, -2.04963448]])



In [7]:

    
plt.contour(corrells)









    Out[7]:





<matplotlib.contour.QuadContourSet at 0x1108da510>



In [8]:

    
dzbin=zdthbin=np.arange(0.002,0.022,0.002)



In [9]:

    
plt.contour(dzbin,zdthbin,corrells)









    Out[9]:





<matplotlib.contour.QuadContourSet at 0x110980e50>



In [10]:

    
ztotal = np.concatenate([np.fliplr(corrells),corrells],axis=1)
ztotal = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal = np.concatenate([-dzbin[::-1],zdthbin],axis=0)
ytotal = np.concatenate([-zdthbin[::-1],zdthbin],axis=0)

plt.figure()
plt.contour(xtotal,ytotal,ztotal)
plt.savefig("../plots/anisotropic2pcfdzzdth.pdf")
plt.show()



In [11]:

    
dr72dat=ascii.read("../output/DR72LCsrarf.dat")



In [12]:

    
dr72dat









    Out[12]:




<Table length=105830>

z ra dec s rar decr
float64 float64 float64 float64 float64 float64
0.450184 56.160452 -0.22136 0.37169 0.980185 -0.003863
0.367165 58.247845 0.216329 0.312739 1.016617 0.003776
0.413324 54.445381 0.619974 0.345944 0.950251 0.010821
0.32247 54.487983 0.485826 0.279501 0.950995 0.008479
0.322646 54.493299 0.48325 0.279634 0.951087 0.008434
0.233469 55.555886 0.582442 0.209831 0.969633 0.010166
0.313645 55.568615 0.608277 0.272806 0.969855 0.010616
0.351716 55.5972 0.434026 0.301375 0.970354 0.007575
0.403208 55.761727 0.495315 0.338761 0.973226 0.008645
0.37691 55.812512 0.521042 0.319842 0.974112 0.009094
... ... ... ... ... ...
0.221777 162.115407 9.12498 0.200306 2.829448 0.159261
0.228503 162.204838 9.136055 0.205796 2.831008 0.159454
0.223359 162.407014 9.242236 0.2016 2.834537 0.161307
0.350793 162.631597 9.112249 0.300692 2.838457 0.159039
0.221194 162.792858 9.065858 0.199829 2.841271 0.158229
0.313262 162.752717 9.235261 0.272514 2.840571 0.161186
0.245009 162.013078 9.587035 0.219143 2.827662 0.167325
0.336967 162.14235 9.941639 0.290404 2.829918 0.173514
0.368113 162.160915 9.861281 0.313432 2.830242 0.172112
0.22131 162.249397 9.950054 0.199924 2.831786 0.173661



In [13]:

    
z=np.mean(dr72dat['z'])



In [14]:

    
z









    Out[14]:





0.32435559778890671



In [18]:

    
om=0.3
ol=0.7
xLCDM=1.0/np.sqrt(om*(1.0+z)**3+ol)
print xLCDM









    



0.846108709066



In [19]:

    
yLCDM0=np.sqrt(om*(1.0+z)**3+ol)/z
print yLCDM0









    



3.64378242875



In [20]:

    
LCDMf = lambda x: 1.0/math.sqrt(om*(1+x)**3+ol)
np.vectorize(LCDMf)

def LCDMfint(z):
    return integrate.quad(LCDMf, 0, z)

LCDMfint=np.vectorize(LCDMfint)

yLCDM=yLCDM0*LCDMfint(z)[0]
print yLCDM









    



1.09180481695



In [21]:

    
xlc=1.0/(1.0+z)

ylc=(1.0+z)/z*np.log(1.0+z)


def ymilne(z):
    return 1.0+z/2.0



In [22]:

    
print xlc
print ylc
print ymilne(z)









    



0.755084209762
1.14703098388
1.16217779889



In [23]:

    
sparLCDM=xLCDM*dzbin
sperLCDM=xLCDM*yLCDM*zdthbin
print sparLCDM
print sperLCDM









    



[ 0.00169222  0.00338443  0.00507665  0.00676887  0.00846109  0.0101533
  0.01184552  0.01353774  0.01522996  0.01692217]
[ 0.00184757  0.00369514  0.00554271  0.00739028  0.00923786  0.01108543
  0.012933    0.01478057  0.01662814  0.01847571]



In [24]:

    
ztotal = np.concatenate([np.fliplr(corrells),corrells],axis=1)
ztotal = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal = np.concatenate([-sparLCDM[::-1],sparLCDM],axis=0)
ytotal = np.concatenate([-sperLCDM[::-1],sperLCDM],axis=0)

plt.figure()
plt.contour(xtotal,ytotal,ztotal)
plt.savefig("../plots/ani2pcfLCDM.pdf")
plt.show()



In [ ]:



In [26]:

    
sparLC=xlc*dzbin
sperLC=xlc*ylc*zdthbin
print sparLC
print sperLC









    



[ 0.00151017  0.00302034  0.00453051  0.00604067  0.00755084  0.00906101
  0.01057118  0.01208135  0.01359152  0.01510168]
[ 0.00173221  0.00346442  0.00519663  0.00692884  0.00866105  0.01039326
  0.01212547  0.01385768  0.01558989  0.0173221 ]



In [27]:

    
ztotal = np.concatenate([np.fliplr(corrells),corrells],axis=1)
ztotal = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal = np.concatenate([-sparLC[::-1],sparLC],axis=0)
ytotal = np.concatenate([-sperLC[::-1],sperLC],axis=0)

plt.figure()
plt.contour(xtotal,ytotal,ztotal)
plt.savefig("../plots/ani2pcfLC.pdf")
plt.show()



In [ ]:



In [28]:

    
sparMilne=xlc*dzbin
sperMilne=xlc*ymilne(z)*zdthbin
print sparMilne
print sperMilne









    



[ 0.00151017  0.00302034  0.00453051  0.00604067  0.00755084  0.00906101
  0.01057118  0.01208135  0.01359152  0.01510168]
[ 0.00175508  0.00351017  0.00526525  0.00702034  0.00877542  0.01053051
  0.01228559  0.01404067  0.01579576  0.01755084]



In [29]:

    
ztotal = np.concatenate([np.fliplr(corrells),corrells],axis=1)
ztotal = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal = np.concatenate([-sparMilne[::-1],sparMilne],axis=0)
ytotal = np.concatenate([-sperMilne[::-1],sperMilne],axis=0)

plt.figure()
plt.contour(xtotal,ytotal,ztotal)
plt.savefig("../plots/ani2pcfMilne.pdf")
plt.show()



In [31]:

    
#ztotal1 = np.concatenate([np.fliplr(corrells),corrells],axis=1)
#ztotal1 = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal1 = np.concatenate([-sparMilne[::-1],sparMilne],axis=0)
ytotal1 = np.concatenate([-sperMilne[::-1],sperMilne],axis=0)
#ztotal2 = np.concatenate([np.fliplr(corrells),corrells],axis=1)
#ztotal2 = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal2 = np.concatenate([-sparLC[::-1],sparLC],axis=0)
ytotal2 = np.concatenate([-sperLC[::-1],sperLC],axis=0)
#ztotal3 = np.concatenate([np.fliplr(corrells),corrells],axis=1)
#ztotal3 = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal3 = np.concatenate([-sparLCDM[::-1],sparLCDM],axis=0)
ytotal3 = np.concatenate([-sperLCDM[::-1],sperLCDM],axis=0)
plt.figure()
plt.contour(xtotal1,ytotal1,ztotal)
plt.contour(xtotal2,ytotal2,ztotal)
plt.contour(xtotal3,ytotal3,ztotal)
plt.savefig("../plots/ani2pcfcompare.pdf")
plt.show()



In [ ]:



In [ ]:



In [ ]:



In [8]:

    
#z=np.arange(0.01,5,0.01)

#om=0.3
#ol=0.7
#yLCDM0=np.sqrt(om*(1+z)**3+ol)/z
#
#LCDMf = lambda x: 1.0/math.sqrt(om*(1+x)**3+ol)
#np.vectorize(LCDMf)
#
#def LCDMfint(zt):
#    return integrate.quad(LCDMf, 0, zt)
#
#LCDMfint=np.vectorize(LCDMfint)
#
#yLCDM=yLCDM0*LCDMfint(z)[0]

#Ez = lambda z: 1/np.sqrt(0.3*(1+z)**3+0.7)

#np.vectorize(Ez)

#def ylcdm(z):
#    for i in range(len(z)):    
#        return 1/z[i]*np.sqrt(0.3*(1+z[i])**3+0.7)/integrate.quad(Ez,0,z[i])
        
#print ylcdm(z)

#ylcdm=1/z*np.

plt.figure()
plt.xlabel('redshift(z)')
plt.ylabel('y(z)')
plt.plot(z,yLCDM,'red',label='yLCDM')
plt.plot(z,ylc,'green',label='yLC')
#plt.plot(z,ylcdm)
plt.plot(z,ymilne(z),'blue',label='yMilne')
plt.legend()
plt.savefig("yz_plots.pdf")
#plt.axes()
#plt.plot(z,ylcdm(z))
plt.show()

200k plots



In [78]:

    
dd2d=np.array([[ 131072.,   37018.,   39251.,   43003.,   46928.,   52380.,
          56755.,   63357.,   70037.,   76252.],
       [  13635.,   25181.,   33553.,   39651.,   45318.,   50737.,
          57137.,   63280.,   69718.,   76424.],
       [   7808.,   18437.,   28142.,   35773.,   43037.,   49519.,
          55848.,   62215.,   69007.,   76304.],
       [   5782.,   15833.,   24651.,   33141.,   40864.,   47585.,
          55253.,   61882.,   68718.,   75510.],
       [   5041.,   13955.,   23161.,   31557.,   38796.,   46402.,
          53552.,   60708.,   67437.,   74405.],
       [   4584.,   13421.,   21805.,   30038.,   37630.,   45206.,
          52575.,   60110.,   66250.,   73244.],
       [   4433.,   13015.,   21007.,   29050.,   36717.,   44679.,
          51354.,   58511.,   65240.,   72592.],
       [   4287.,   12526.,   20387.,   28567.,   36093.,   43421.,
          50931.,   57911.,   64958.,   71831.],
       [   4093.,   12155.,   20054.,   28384.,   35523.,   43465.,
          50161.,   56795.,   64062.,   70675.],
       [   4170.,   11915.,   19636.,   27596.,   35199.,   42741.,
          49570.,   56936.,   63389.,   70143.]])
rr2d=np.array([[ 219389.,   52749.,   86685.,  118306.,  150065.,  178945.,
         209520.,  239222.,  267225.,  297534.],
       [  17753.,   52569.,   86627.,  117295.,  148557.,  178415.,
         208245.,  238385.,  267604.,  294867.],
       [  17566.,   51819.,   85165.,  116994.,  147427.,  177454.,
         207106.,  236741.,  264963.,  293735.],
       [  17780.,   52130.,   84770.,  115177.,  146679.,  176461.,
         205271.,  233772.,  263140.,  291924.],
       [  17734.,   51545.,   83932.,  114524.,  145841.,  175503.,
         204296.,  232320.,  260856.,  289001.],
       [  17313.,   51400.,   83842.,  114798.,  145069.,  173206.,
         203071.,  231689.,  259312.,  287118.],
       [  17095.,   50846.,   82932.,  113880.,  143613.,  173206.,
         201398.,  228866.,  258066.,  286164.],
       [  16866.,   50369.,   82171.,  112970.,  142838.,  171366.,
         199218.,  227239.,  254381.,  283645.],
       [  16956.,   49662.,   81561.,  112150.,  141164.,  170538.,
         198573.,  226388.,  253592.,  280492.],
       [  16865.,   49297.,   81354.,  111885.,  140364.,  168967.,
         196672.,  223977.,  251690.,  278476.]])
dr2d=np.array([[  17225.,   51016.,   83901.,  115532.,  145602.,  174413.,
         205404.,  235251.,  261629.,  291459.],
       [  17130.,   50742.,   83297.,  114374.,  145734.,  174830.,
         204178.,  232755.,  260633.,  289584.],
       [  16806.,   50308.,   82561.,  114032.,  144563.,  173730.,
         203012.,  231303.,  260188.,  288244.],
       [  16964.,   50424.,   82672.,  113163.,  143783.,  173336.,
         202119.,  231121.,  258214.,  285733.],
       [  16846.,   49890.,   81972.,  113003.,  143170.,  172195.,
         201662.,  229361.,  256191.,  284896.],
       [  16966.,   49739.,   81756.,  112037.,  141410.,  170703.,
         199887.,  227671.,  256364.,  283776.],
       [  16578.,   49172.,   80573.,  111373.,  141073.,  170247.,
         198875.,  225671.,  254128.,  281677.],
       [  16713.,   48801.,   80484.,  110405.,  139692.,  169457.,
         197079.,  225067.,  252907.,  280422.],
       [  16613.,   49008.,   80096.,  109819.,  139137.,  167602.,
         194844.,  223229.,  251104.,  278213.],
       [  16212.,   47944.,   79251.,  109194.,  138363.,  165319.,
         194678.,  222679.,  248924.,  276076.]])



In [79]:

    
correl200k=(dd2d-4.0*dr2d+4.0*rr2d)/rr2d



In [80]:

    
dzbin









    Out[80]:





array([ 0.002,  0.004,  0.006,  0.008,  0.01 ,  0.012,  0.014,  0.016,
        0.018,  0.02 ])



In [82]:

    
ztotal = np.concatenate([np.fliplr(correl200k),correl200k],axis=1)
ztotal = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal = np.concatenate([-dzbin[::-1],zdthbin],axis=0)
ytotal = np.concatenate([-zdthbin[::-1],zdthbin],axis=0)

plt.figure()
plt.contour(xtotal,ytotal,ztotal)
plt.savefig("../plots/anisotropic2pcf200k.pdf")
plt.show()



In [83]:

    
correl200k









    Out[83]:





array([[ 4.28338704,  0.83319115,  0.5812655 ,  0.45728027,  0.43167961,
         0.39402051,  0.34946067,  0.33124462,  0.34585462,  0.33795129],
       [ 0.90840985,  0.61802583,  0.54108996,  0.43765719,  0.38106585,
         0.36475072,  0.35249346,  0.35992197,  0.36472549,  0.33084747],
       [ 0.61755664,  0.47243289,  0.45274467,  0.40703797,  0.369627  ,
         0.36299548,  0.34872964,  0.35467874,  0.33252567,  0.33454644],
       [ 0.5087739 ,  0.43462498,  0.38979592,  0.35768426,  0.35756993,
         0.34050017,  0.33059224,  0.31007135,  0.33602645,  0.34349351],
       [ 0.48454945,  0.39916578,  0.36935853,  0.32867347,  0.33927359,
         0.33978906,  0.31370169,  0.31225895,  0.33005566,  0.31427227],
       [ 0.34494311,  0.39036965,  0.35959304,  0.35786338,  0.36028373,
         0.31879958,  0.32161658,  0.32881147,  0.30095792,  0.30165994],
       [ 0.38028663,  0.38766078,  0.36708388,  0.34315068,  0.32641195,
         0.32628777,  0.30509737,  0.31149668,  0.3138422 ,  0.316392  ],
       [ 0.29046603,  0.37320574,  0.33022599,  0.34369302,  0.34078467,
         0.29794125,  0.29860254,  0.29307909,  0.27853495,  0.29869379],
       [ 0.32230479,  0.29743063,  0.31772538,  0.33622827,  0.30908022,
         0.3237343 ,  0.32772331,  0.30669028,  0.29186252,  0.28446801],
       [ 0.4021346 ,  0.35148183,  0.34476485,  0.34285204,  0.3077926 ,
         0.33931478,  0.29259884,  0.27738562,  0.29581231,  0.28635502]])



In [85]:

    
ztotal = np.concatenate([np.fliplr(correl200k),correl200k],axis=1)
ztotal = np.concatenate([np.flipud(ztotal),ztotal],axis=0)
xtotal = np.concatenate([-dzbin[::-1],zdthbin],axis=0)
ytotal = np.concatenate([-zdthbin[::-1],zdthbin],axis=0)

plt.figure()
plt.contour(xtotal,ytotal,ztotal,levels=(0.2,0.3,0.4,0.5,0.6,0.7,0.8))
plt.savefig("../plots/anisotropic2pcf200kl.pdf")
plt.show()



In [ ]:



In [64]:

    
dzmu1









    Out[64]:





array([ 0.002,  0.004,  0.006,  0.008,  0.01 ,  0.012,  0.014,  0.016,
        0.018,  0.02 ])



In [65]:

    
correlmu1









    Out[65]:





array([ 2.08273679, -0.06591932, -1.14525488, -1.49923222, -1.74361248,
       -1.81418446, -1.90115591, -1.94329043, -1.93130228, -1.96064259])



In [100]:

    
sLCDM=pairdats['xLCDM']*np.sqrt((pairdats['deltaZ'])**2+(pairdats['yLCDM']*pairdats['Zdeltatheta'])**2)

z	ra	dec	s	rar	decr
float64	float64	float64	float64	float64	float64
0.450184	56.160452	-0.22136	0.37169	0.980185	-0.003863
0.367165	58.247845	0.216329	0.312739	1.016617	0.003776
0.413324	54.445381	0.619974	0.345944	0.950251	0.010821
0.32247	54.487983	0.485826	0.279501	0.950995	0.008479
0.322646	54.493299	0.48325	0.279634	0.951087	0.008434
0.233469	55.555886	0.582442	0.209831	0.969633	0.010166
0.313645	55.568615	0.608277	0.272806	0.969855	0.010616
0.351716	55.5972	0.434026	0.301375	0.970354	0.007575
0.403208	55.761727	0.495315	0.338761	0.973226	0.008645
0.37691	55.812512	0.521042	0.319842	0.974112	0.009094
...	...	...	...	...	...
0.221777	162.115407	9.12498	0.200306	2.829448	0.159261
0.228503	162.204838	9.136055	0.205796	2.831008	0.159454
0.223359	162.407014	9.242236	0.2016	2.834537	0.161307
0.350793	162.631597	9.112249	0.300692	2.838457	0.159039
0.221194	162.792858	9.065858	0.199829	2.841271	0.158229
0.313262	162.752717	9.235261	0.272514	2.840571	0.161186
0.245009	162.013078	9.587035	0.219143	2.827662	0.167325
0.336967	162.14235	9.941639	0.290404	2.829918	0.173514
0.368113	162.160915	9.861281	0.313432	2.830242	0.172112
0.22131	162.249397	9.950054	0.199924	2.831786	0.173661