In [1]:
import healpix_util as hu
import astropy as ap
import math, sys, time
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
from astropy.stats import bayesian_blocks
from astropy.stats import histogram
from astropy.visualization import hist
%matplotlib inline
#%precision 4
#plt.style.use('ggplot')
In [2]:
hmqdat1=ascii.read("./hmqdata_1.dat")
Z=sum(hmqdat1['Z'])/len(hmqdat1['Z'])
print Z
In [3]:
NSIDE=512
hmqhp=hu.HealPix("RING",NSIDE)
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)
newcol=hmqdat1.Column(name='PIX',data=hmqhppix)
hmqdat1.add_column(newcol)
In [4]:
#hmqdata1p = open("hmqdata12_fp1.dat",'w')
#hmqdata12p.write(hmqdat12)
#ascii.write(hmqdat1,output="hmqdata12_fp1.dat")
hpixangth,hpixangphi=hu.pix2ang(NSIDE,hmqdat1[:]['PIX'])
In [ ]:
dddata = open("./output/DD_1ff_simpler.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
Z1=hmqdat1[0]['Z']
while j<len(hmqdat1['PIX']):
try:
Z2=hmqdat1[j]['Z']
deltaZ=abs(Z1-Z2)
if deltaZ<0.02:
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
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)
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))
else:
break
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()
In [10]:
abs(0.5)
Out[10]:
In [6]:
hmqdat1
Out[6]:
In [ ]: