<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Healpix pixelization of DR72 SDSS Database
First import all the modules such as healpy and astropy needed for analyzing the structure
Step1:  Read the data file
Sorted and reduced column set data can now be 'read' to reduce RAM requirements of the table reading.
Step2:  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.
Step3:  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.
Step4:  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.
<ASSISTANT_TASK:>
Python Code:
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
Explanation: Healpix pixelization of DR72 SDSS Database
First import all the modules such as healpy and astropy needed for analyzing the structure
End of explanation
sdssdr72=ascii.read('/home/rohin/Desktop/healpix/sdssdr72_sorted_z.dat')
Explanation: Read the data file
Sorted and reduced column set data can now be 'read' to reduce RAM requirements of the table reading.
End of explanation
NSIDE=64
dt72hpix=hu.HealPix("ring",NSIDE)
Explanation: 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.
End of explanation
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()
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()
Explanation: 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.
End of explanation
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
hpixdata
hu.orthview(hpixdata,rot=180)
pixcl=hu.anafast(hpixdata,lmax=300)
ell = np.arange(len(pixcl))
plt.figure()
plt.plot(ell,np.log(pixcl))
plt.show()
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()
theta=np.arange(0,np.pi,0.001)
correldat = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(pixcl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correldat[0:600])
plt.show()
plt.figure()
plt.plot(theta*180/math.pi,correldat)
plt.show()
randra,randdec=hu.randsphere(2200000)
randhp=hu.HealPix("RING",NSIDE)
randhppix=randhp.eq2pix(randra,randdec)
randpixdat=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(randhppix)):
    randpixdat[randhppix[j]]+=1
randmaphp=hu.mollview(randpixdat)
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()
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()
finalcorrel=correldat-correlrand
plt.figure()
plt.plot(theta[0:600]*180/math.pi,finalcorrel[0:600])
plt.show()
finalpix=hpixdata-randpixdat
hu.mollview(finalpix,rot=180)
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()
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()
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()
Explanation: 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.
End of explanation
<END_TASK>