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
from scipy import integrate


    

In [16]:

    

hmqdat12=ascii.read("/home/rohin/Documents/ipy_notebooks/healpix_practice/hmqdata_12.dat")

Z=sum(hmqdat12['Z'])/len(hmqdat12['Z'])
print Z


    

    




5.85858536585

In [17]:

    

NSIDE=128
hmqhp=hu.HealPix("RING",NSIDE)


    

In [18]:

    

hmqhppix=hmqhp.eq2pix(hmqdat12['RA'],hmqdat12['DEC'])
hmqpixdat=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(hmqhppix)):
    hmqpixdat[hmqhppix[j]]+=1
hu.orthview(hmqpixdat)


    

In [19]:

    

newcol=hmqdat12.Column(name='PIX',data=hmqhppix)
hmqdat12.add_column(newcol)


    

In [20]:

    

hmqdata12p = open("hmqdata1_fp12.dat",'w')
#hmqdata12p.write(hmqdat12)
ascii.write(hmqdat12,output="hmqdata12_fp1.dat")
hpixangth,hpixangphi=hu.pix2ang(NSIDE,hmqdat12[:]['PIX'])


    

In [21]:

    

dddata = open("./DD_12.dat",'w')
dddata.write("Z1\t RA1\t DEC1\t Z2\t RA2\t DEC2\t deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\t thetadd \t DD \n")

while len(hmqdat12['PIX'])>0:
    j=1
    while j<len(hmqdat12['PIX']):
        try:
            Z1=hmqdat12[0]['Z']
            Z2=hmqdat12[j]['Z']
            RA1=hmqdat12[0]['RA']*math.pi/180.0
            RA2=hmqdat12[j]['RA']*math.pi/180.0
            DEC1=hmqdat12[0]['DEC']*math.pi/180.0
            DEC2=hmqdat12[j]['DEC']*math.pi/180.0
            deltaZ=Z1-Z2
            dddata.write("%f\t" %Z1)
            dddata.write("%f\t" %RA1)
            dddata.write("%f\t" %DEC1)
            dddata.write("%f\t" %Z2)
            dddata.write("%f\t" %RA2)
            dddata.write("%f\t" %DEC2)
            deltatheta=math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2))
            Zdeltatheta=Z*deltatheta
            R=math.sqrt(deltaZ**2+Zdeltatheta**2)
            alpha=math.acos(deltaZ/R)
            dddata.write("%f\t" %deltaZ)
            dddata.write("%f\t" %deltatheta)
            dddata.write("%f\t" %Zdeltatheta)
            dddata.write("%f\t" %R)
            dddata.write("%f\t" %alpha)
            DD=hmqpixdat[hmqdat12[0]['PIX']]*hmqpixdat[hmqdat12[j]['PIX']]
            thetadd=math.acos(math.cos(hpixangth[0])*math.cos(hpixangth[j])+math.sin(hpixangth[0])*math.sin(hpixangth[j])*math.cos(hpixangphi[0]-hpixangphi[j]))
            dddata.write("%f \t %f \n"%(thetadd,DD))
        except ValueError as e:
            print hpixangth[0]
            print hpixangth[j]
            print hpixangphi[0]
            print hpixangphi[j]
        j=j+1
    hmqdat12.remove_row(0)
dddata.close()


    

In [12]:

    

dddata = open("./DD_12f.dat",'w')
dddata.write("Z1\t RA1\t DEC1\t Z2\t RA2\t DEC2\t deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\t thetadd \t DD \n")

while len(hmqdat12['PIX'])>0:
    j=1
    while j<len(hmqdat12['PIX']):
        try:
            Z1=hmqdat12[0]['Z']
            Z2=hmqdat12[j]['Z']
            RA1=hmqdat12[0]['RA']*math.pi/180.0
            RA2=hmqdat12[j]['RA']*math.pi/180.0
            DEC1=hmqdat12[0]['DEC']*math.pi/180.0
            DEC2=hmqdat12[j]['DEC']*math.pi/180.0
            deltaZ=Z1-Z2
            deltatheta=math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2))
            Zdeltatheta=Z*deltatheta
            R=math.sqrt(deltaZ**2+Zdeltatheta**2)
            alpha=math.acos(deltaZ/R)
            if deltaZ<0.2:
                dddata.write("%f\t" %Z1)
                dddata.write("%f\t" %RA1)
                dddata.write("%f\t" %DEC1)
                dddata.write("%f\t" %Z2)
                dddata.write("%f\t" %RA2)
                dddata.write("%f\t" %DEC2)
                dddata.write("%f\t" %deltaZ)
                dddata.write("%f\t" %deltatheta)
                dddata.write("%f\t" %Zdeltatheta)
                dddata.write("%f\t" %R)
                dddata.write("%f\t" %alpha)
                DD=hmqpixdat[hmqdat12[0]['PIX']]*hmqpixdat[hmqdat12[j]['PIX']]
                thetadd=math.acos(math.cos(hpixangth[0])*math.cos(hpixangth[j])+math.sin(hpixangth[0])*math.sin(hpixangth[j])*math.cos(hpixangphi[0]-hpixangphi[j]))
                dddata.write("%f \t %f \n"%(thetadd,DD))
        except ValueError as e:
            print hpixangth[0]
            print hpixangth[j]
            print hpixangphi[0]
            print hpixangphi[j]
        j=j+1
    hmqdat12.remove_row(0)
dddata.close()


    

In [23]:

    

#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(dddata[:]['R'],dddata[:]['DD'])
plt.xlabel('R')
plt.ylabel('DD')
plt.show()


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-23-4b51179a32d5> in <module>()
      2 #correldat = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(pixcl))/(4*math.pi)
      3 plt.figure()
----> 4 plt.plot(dddata[:]['R'],dddata[:]['DD'])
      5 plt.xlabel('R')
      6 plt.ylabel('DD')

TypeError: 'file' object has no attribute '__getitem__'





    






<matplotlib.figure.Figure at 0x7f15eee84d10>

In [27]:

    

hmqdat12['PIX']


    

    
Out[27]:




<Column name='PIX' dtype='int64' length=0>

In [ ]: