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
import pp
import gnumpy as gnp
from scipy import integrate
%matplotlib inline
#%precision 4
#plt.style.use('ggplot')


    

In [3]:

    

hmqdat1=ascii.read("./hmqdata_1.dat")

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


    

    




0.238083091217

In [4]:

    

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


    

In [5]:

    

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


    

In [6]:

    

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


    

In [7]:

    

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


    

    




WARNING: AstropyDeprecationWarning: hmqdata1_fp1.dat already exists. Automatically overwriting ASCII files is deprecated. Use the argument 'overwrite=True' in the future. [astropy.io.ascii.ui]

In [14]:

    

dddata = open("./DD_1ff.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(hmqdat1['PIX'])>0:
    j=1
    while j<len(hmqdat1['PIX']):
        try:
            Z1=hmqdat1[0]['Z']
            Z2=hmqdat1[j]['Z']
            RA1=hmqdat1[0]['RA']*math.pi/180.0
            RA2=hmqdat1[j]['RA']*math.pi/180.0
            DEC1=hmqdat1[0]['DEC']*math.pi/180.0
            DEC2=hmqdat1[j]['DEC']*math.pi/180.0
            deltaZ=abs(Z1-Z2)
            deltatheta=abs(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 and Zdeltatheta<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[hmqdat1[0]['PIX']]*hmqpixdat[hmqdat1[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
    hmqdat1.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 [10]:

    

abs(0.5)


    

    
Out[10]:





0.5

In [ ]: