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 [2]:

    

hmqdat1=ascii.read("/home/rohin/Documents/ipy_notebooks/healpix_practice/hmqdata_1.dat")

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


    

    




0.238083091217

In [3]:

    

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


    

In [4]:

    

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)


    

    




/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:190: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  a = empty(shape, dtype, order)






    

In [5]:

    

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


    

In [6]:

    

hmqdata1p = open("hmqdata1_fp1.dat",'w')
#hmqdata12p.write(hmqdat12)
ascii.write(hmqdat1,output="hmqdata1_fp1.dat")


    

In [7]:

    

hu.mollview(hmqpixdat,rot=180,norm='hist')


    

    




/home/rohin/.local/lib/python2.7/site-packages/healpy-1.8.6-py2.7-linux-x86_64.egg/healpy/projaxes.py:808: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.
Check the NumPy 1.11 release notes for more information.
  result[np.isinf(val.data)] = -np.inf






    

In [10]:

    

hpixangth,hpixangphi=hu.pix2ang(NSIDE,hmqdat1[:]['PIX'])


    

In [ ]:

    

dddata = open("./DD_1.dat",'w')
dddata.write("thetadd \t DD \n")
#while len(hmqdat1['PIX']))>0:
while len(hmqdat1['PIX'])>0:
    j=1
    while j<len(hmqdat1['PIX']):
        try:
            DD=hmqpixdat[hmqdat1[1]['PIX']]*hmqpixdat[hmqdat1[j]['PIX']]
            thetadd=math.acos(math.cos(hpixangth[1])*math.cos(hpixangth[j])+math.sin(hpixangth[1])*math.sin(hpixangth[j])*math.cos(hpixangphi[1]-hpixangphi[j]))
            dddata.write("%f \t %f \n"%(thetadd,DD))
        except ValueError as e:
            print hpixangth[1]
            print hpixangth[j]
            print hpixangphi[1]
            print hpixangphi[j]
        j=j+1
    hmqdat1.remove_row(0)
dddata.close()


    

In [9]:

    

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.orthview(randpixdat)
hpixangth,hpixangphi=hu.pix2ang(NSIDE,randpixdat[:])


    

In [ ]:

    

nndata = open("./NN_1.dat",'w')
nndata.write("thetann \t NN \n")
for i in range(len(hpixdata)):
    for j in range(len(hpixdata)):
        try:
            NN=hpixdata[pixdata[i]['pix']]*hpixdata[pixdata[j]['pix']]
            #theta,phi=hu.pix2ang(NSIDE,pixdata
            thetann=math.acos(math.cos(hpixangth[i])*math.cos(hpixangth[j])+math.sin(hpixangth[i])*math.sin(hpixangth[j])*math.cos(hpixangphi[i]-hpixangphi[j]))
            nndata.write("%f \t %f \n"%(thetann,NN))
        except ValueError as e:
            print hpixangth[i]
            print hpixangth[j]
            print hpixangphi[i]
            print hpixangphi[j]
nndata.close()


    

In [12]:

    

hmqdat1['PIX']


    

    
Out[12]:




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

189724
573232
1243348
1393677
2315306
2202681
2382947
2385041
1021181
2344192
1740039
523639
...
957399
1663962
1227741
1469405
431492
1549282
1625060
1819621
1641454
1137647
1547248
1770489

In [2]:

    

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 [9]:

    

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


    

In [10]:

    

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 [7]:

    

sum(hmqpixdat)


    

    
Out[7]:





41.0

In [8]:

    

hmqpixdat


    

    
Out[8]:





array([ 0.,  0.,  0., ...,  0.,  0.,  0.])

In [16]:

    

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


    

In [19]:

    

dddata = open("./DD_12.dat",'w')
dddata.write("thetadd \t DD \n")
#while len(hmqdat1['PIX']))>0:
while len(hmqdat12['PIX'])>0:
    j=1
    while j<len(hmqdat12['PIX']):
        try:
            DD=hmqpixdat[hmqdat12[1]['PIX']]*hmqpixdat[hmqdat12[j]['PIX']]
            thetadd=math.acos(math.cos(hpixangth[1])*math.cos(hpixangth[j])+math.sin(hpixangth[1])*math.sin(hpixangth[j])*math.cos(hpixangphi[1]-hpixangphi[j]))
            dddata.write("%f \t %f \n"%(thetadd,DD))
        except ValueError as e:
            print hpixangth[1]
            print hpixangth[j]
            print hpixangphi[1]
            print hpixangphi[j]
        j=j+1
    hmqdat12.remove_row(0)
dddata.close()


    

In [ ]: