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.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 [4]:
sdssq=ascii.read("/home/rohin/Desktop/healpix/SDSS_quasar_radecz.dat")
In [5]:
sdssq
Out[5]:
In [6]:
sdssq.sort("z")
In [40]:
newcol=sdssq.Column(name='pix',data=pixdata)
sdssq.add_column(newcol)
Out[40]:
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 [19]:
NSIDE=128
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 [20]:
j=0
for i in range(1,25):
pixdata = open("/home/rohin/Desktop/healpix/binned3/pixdata%d_%d.dat"%(NSIDE,i),'w')
pixdata.write("ra\t dec\t z\t pix \n")
#for j in range(len(sdssdr72)):
while sdssq[j]['z']<0.2*i:
pixdata.write("%f\t" %sdssq[j]['RA'])
pixdata.write("%f\t" %sdssq[j]['Dec'])
pixdata.write("%f\t" %sdssq[j]['z'])
pixdata.write("%d\n" %dt72hpix.eq2pix(sdssq[j]['RA'],sdssq[j]['Dec']))
#print dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec'])
j=j+1
pixdata.close()
In [21]:
pixdata = ascii.read("/home/rohin/Desktop/healpix/binned3/pixdata128_9.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [22]:
hpixdata
Out[22]:
In [23]:
hu.mollview(hpixdata,rot=180)
In [25]:
blah=hu.read_map("/home/rohin/Desktop/healpix/rast_sdss_dr72safe0_nside512.fits")
In [27]:
hu.mollview(blah,rot=180)
In [ ]: