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 [66]:
hmqdat12=ascii.read("/home/rohin/Documents/ipy_notebooks/healpix_practice/hmqdata_12.dat")

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

print "Z=%f\t" %Z

hmqdata12a = open("hmqdata12_analysed.dat",'w')
hmqdata12a.write("Z1\t RA1\t DEC1\t Z2\t RA2\t DEC2\t deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\n")

while len(hmqdat12['Z'])>0:
    j=1
    while j<len(hmqdat12['Z']):
        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
        hmqdata12a.write("%f\t" %Z1)
        hmqdata12a.write("%f\t" %RA1)
        hmqdata12a.write("%f\t" %DEC1)
        hmqdata12a.write("%f\t" %Z2)
        hmqdata12a.write("%f\t" %RA2)
        hmqdata12a.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)
        hmqdata12a.write("%f\t" %deltaZ)
        hmqdata12a.write("%f\t" %deltatheta)
        hmqdata12a.write("%f\t" %Zdeltatheta)
        hmqdata12a.write("%f\t" %R)
        hmqdata12a.write("%f\n" %alpha)
        j=j+1
    hmqdat12.remove_row(0)
hmqdata12a.close()


Z=5.858585	

In [64]:
theta1,phi1=hu.eq2ang(RA1,DEC1)
print RA1
print DEC1
print theta1
print phi1


340.204083
-18.662167
1.8965125865
5.93768137708

In [58]:
help(hu.eq2ang)


Help on function eq2ang in module healpix_util.coords:

eq2ang(ra, dec)
    convert equatorial ra,dec in degrees to angular theta, phi in radians
    
    parameters
    ----------
    ra: scalar or array
        Right ascension in degrees
    dec: scalar or array
        Declination in degrees
    
    returns
    -------
    theta,phi: tuple
        theta = pi/2-dec*D2R # in [0,pi]
        phi   = ra*D2R       # in [0,2*pi]


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

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

print "Z=%f\t" %Z

hmqdata1a = open("hmqdata1_analysed.dat",'w')
hmqdata1a.write("Z1\t RA1\t DEC1\t Z2\t RA2\t DEC2\t deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\n")

while len(hmqdat1['Z'])>0:
    j=1
    while j<len(hmqdat1['Z']):
        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=Z1-Z2
        hmqdata1a.write("%f\t" %Z1)
        hmqdata1a.write("%f\t" %RA1)
        hmqdata1a.write("%f\t" %DEC1)
        hmqdata1a.write("%f\t" %Z2)
        hmqdata1a.write("%f\t" %RA2)
        hmqdata1a.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)
        hmqdata1a.write("%f\t" %deltaZ)
        hmqdata1a.write("%f\t" %deltatheta)
        hmqdata1a.write("%f\t" %Zdeltatheta)
        hmqdata1a.write("%f\t" %R)
        hmqdata1a.write("%f\n" %alpha)
        j=j+1
    hmqdat1.remove_row(0)
hmqdata1a.close()

In [ ]: