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
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 cython_metric import *
from progressbar import *
#from astroML.datasets import fetch_sdss_specgals
#from astroML.correlation import bootstrap_two_point_angular
%matplotlib inline
Random catalog created based on the survey limitations also taken from http://cosmo.nyu.edu/~eak306/SDSS-LRG.html
In [2]:
dataR=fits.open("/Users/rohin/Downloads/random-DR7-Full.fits")
In [3]:
dataR
Out[3]:
In [4]:
dataR=dataR[1].data
In [5]:
len(dataR)
Out[5]:
In [6]:
dataR
Out[6]:
In [7]:
dataR.columns
Out[7]:
In [8]:
dataR['Z']
Out[8]:
In [9]:
dataR['RA']
Out[9]:
In [10]:
dataR['DEC']
Out[10]:
In [11]:
Ez = lambda x: 1/m.sqrt(0.3*(1+x)**3+0.7)
np.vectorize(Ez)
#Calculate comoving distance of a data point using the Redshift - This definition is based on the cosmology model we take. Here the distance for E-dS universe is considered. Also note that c/H0 ratio is cancelled in the equations and hence not taken.
def DC_LCDM(z):
return integrate.quad(Ez, 0, z)[0]
DC_LCDM=np.vectorize(DC_LCDM)
In [14]:
rdr7f = open("./output/randDR7srarf.dat",'w')
rdr7f.write("z\t ra\t dec\t s\t rar\t decr \n")
for i in range(0,len(dataR)-1):
rdr7f.write("%f\t " %dataR['Z'][i])
rdr7f.write("%f\t %f\t " %(dataR['RA'][i],dataR['DEC'][i]))
rdr7f.write("%f\t " %DC_LCDM(dataR['Z'][i]))
rdr7f.write("%f\t %f\n " %(dataR['RA'][i]*pi/180.0,dataR['DEC'][i]*pi/180.0))
rdr7f.close()
In [15]:
data=ascii.read("./output/randDR7srarf.dat")
In [16]:
data
Out[16]:
In [17]:
len(data)
Out[17]:
In [18]:
data.remove_column('z')
data.remove_column('ra')
data.remove_column('dec')
In [19]:
data
Out[19]:
In [20]:
s=np.array(data['s'])
rar=np.array(data['rar'])
decr=np.array(data['decr'])
In [21]:
dat=np.array([s,rar,decr])
In [22]:
dat.reshape(3,len(data['s']))
Out[22]:
In [23]:
dat=dat.transpose()
In [24]:
dat
Out[24]:
In [25]:
LCDMmetric(dat[0],dat[1])
Out[25]:
In [26]:
bins=np.arange(0.,0.08,0.005)
In [27]:
print bins
In [28]:
Nbins=len(bins)
In [29]:
Nbins
Out[29]:
In [30]:
BT_R = BallTree(dat,metric='pyfunc',func=LCDMmetric)
temp = pickle.dumps(BT_R)
BTR = pickle.loads(temp)
In [ ]:
pb=ProgressBar()
In [ ]:
start_time=time.time()
@pickle_results("RDR72RR.pkl")
def rrcal(BTR,dat,bins,Nbins):
counts_RR=np.zeros(Nbins)
for i in pb(range(Nbins)):
counts_RR[i]=np.sum(BTR.query_radius(dat, bins[i],count_only=True))
RR = np.diff(counts_RR)
print counts_RR
print RR
return RR
#def mf_wrap(args):
# return ddcal(*args)
#pool=mp.Pool(8)
#arg=[(BTD,dat,bins,Nbins)]
RRresult = rrcal(BTR,dat,bins,Nbins)
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
#DDresult=pool.map(mf_wrap,arg)
In [ ]:
plt.plot(bins[1:len(bins)],RRresult[0],'ro')
In [ ]:
RRresult[0]
In [ ]:
In [ ]: