First import all the modules such as healpy and astropy needed for analyzing the structure
In [172]:
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.constants import c
import matplotlib.pyplot as plt
import math
import scipy.special as sp
Read the data file
Sorted and reduced column set data can now be 'read' to reduce RAM requirements of the table reading.
In [173]:
sdssdr72=ascii.read('/home/rohin/Desktop/healpix/sdssdr72_sorted_z.dat')
Create a healpix map with NSIDE=64 (no. of pixels = 49152 as $NPIX=12\times NSIDE^2$) because the no. of galaxies in the survey are less. For higher resolution (later for dr12) we will consider NSIDE=512 or even 1024. For now, we will create a 64 NSIDE map.
In [174]:
NSIDE=64
dt72hpix=hu.HealPix("ring",NSIDE)
We have data of galaxies with redshifts between 0 and 0.5 ($0<z<0.5$). To look at a time slice/at a certain epoch we need to choose the list of galaxies within a redshift window. As, measurement of redshift has $\pm 0.05$ error. We can bin all the data into redshifts with range limited to 0.05 variation each. So, we have 10 databins with almost identical redshifts. We save each databin in a different file.
In [175]:
j=0
for i in range(1,17):
pixdata = open("/home/rohin/Desktop/healpix/binned1/pixdata%d_%d.dat"%(NSIDE,i),'w')
pixdata.write("ra\t dec\t z\t pix \n")
#for j in range(len(sdssdr72)):
try:
while sdssdr72[j]['z']<0.03*i:
pixdata.write("%f\t" %sdssdr72[j]['ra'])
pixdata.write("%f\t" %sdssdr72[j]['dec'])
pixdata.write("%f\t" %sdssdr72[j]['z'])
pixdata.write("%d\n" %dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec']))
#print dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec'])
j=j+1
except:
pass
pixdata.close()
In [176]:
for i in range(1,17):
pixdata = ascii.read("/home/rohin/Desktop/healpix/binned1/pixdata%d_%d.dat"%(NSIDE,i))
mpixdata = open("/home/rohin/Desktop/healpix/binned1/masked/pixdata%d_%d.dat"%(NSIDE,i),'w')
mpixdata.write("ra\t dec\t z\t pix \n")
for j in range((len(pixdata)-1)):
if 100<pixdata[j]['ra']<250:
mpixdata.write("%f\t" %pixdata[j]['ra'])
mpixdata.write("%f\t" %pixdata[j]['dec'])
mpixdata.write("%f\t" %pixdata[j]['z'])
mpixdata.write("%d\n" %pixdata[j]['pix'])
#pixdata.write("/home/rohin/Desktop/healpix/binned1/masked/pixdata_%d.dat"%i,format='ascii')
#print dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec'])
mpixdata.close()
We now, take each databin and assign the total no. of galaxies as the value of each pixel. The following routine will calculate the no. of galaxies by couting the occurence of pixel numbers in the file.
In [234]:
pixdata = ascii.read("/home/rohin/Desktop/healpix/binned1/masked/pixdata%d_2.dat"%NSIDE)
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [235]:
hpixdata
Out[235]:
In [236]:
hu.orthview(hpixdata,rot=180)
In [213]:
pixcl=hu.anafast(hpixdata,lmax=300)
ell = np.arange(len(pixcl))
plt.figure()
plt.plot(ell,np.log(pixcl))
plt.show()
In [214]:
pixcl=hu.anafast(hpixdata,lmax=300)
ell = np.arange(len(pixcl))
plt.figure()
plt.plot(ell,np.sqrt(ell*(ell+1)*pixcl/(4*math.pi)))
plt.show()
In [215]:
theta=np.arange(0,np.pi,0.001)
In [216]:
correldat = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(pixcl))/(4*math.pi)
In [217]:
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correldat[0:600])
plt.show()
In [218]:
plt.figure()
plt.plot(theta*180/math.pi,correldat)
plt.show()
In [219]:
randra,randdec=hu.randsphere(2200000)
In [220]:
randhp=hu.HealPix("RING",NSIDE)
In [221]:
randhppix=randhp.eq2pix(randra,randdec)
In [222]:
randpixdat=np.array(np.zeros(hu.nside2npix(NSIDE)))
In [223]:
for j in range(len(randhppix)):
randpixdat[randhppix[j]]+=1
In [224]:
randmaphp=hu.mollview(randpixdat)
In [225]:
randcl=hu.anafast(randpixdat,lmax=300)
ell = np.arange(len(randcl))
plt.figure()
plt.plot(ell,np.sqrt(ell*(ell+1)*randcl/(4*math.pi)))
plt.show()
In [226]:
correlrand = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(randcl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correlrand[0:600])
plt.show()
In [227]:
finalcorrel=correldat-correlrand
plt.figure()
plt.plot(theta[0:600]*180/math.pi,finalcorrel[0:600])
plt.show()
In [228]:
finalpix=hpixdata-randpixdat
In [229]:
hu.mollview(finalpix,rot=180)
In [230]:
cl=hu.anafast(finalpix,lmax=300)
ell = np.arange(len(cl))
plt.figure()
plt.plot(ell,np.sqrt(ell*(ell+1)*cl/(4*math.pi)))
plt.show()
In [231]:
correlrand = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(cl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correlrand[0:600])
plt.show()
In [232]:
finalcl=pixcl-randcl
correlrand = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(finalcl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correlrand[0:600])
plt.show()