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
"""


sdssq=ascii.read("/home/rohin/Desktop/healpix/SDSS_quasar_radecz.dat")

sdssq

sdssq.sort("z")

newcol=sdssq.Column(name='pix',data=pixdata)
sdssq.add_column(newcol)

"""
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=128
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,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()

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

hpixdata

hu.mollview(hpixdata,rot=180)

blah=hu.read_map("/home/rohin/Desktop/healpix/rast_sdss_dr72safe0_nside512.fits")

hu.mollview(blah,rot=180)

"""
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
"""


