

In [79]:

    
import numpy as np
import scipy as sp
from scipy import integrate
from math import *
import matplotlib.pyplot as plt
from astropy.io import ascii
import healpix_util as hu

Define angular diameter distances for Flat LambdaCDM cosmology with $\Omega_M=0.3$ and $\Omega_{\Lambda}=0.7$ starting with definition of $E(z)=\sqrt{0.3(1+z)^3+0.7}$. As we know Angular diameter distance ($D_A$) is comoving distance ($D_C$) divided by $1+z$. We write $D_A=\frac{D_C}{1+z}=\frac{1}{1+z}\int_0^z\frac{dz}{E(z)}$



In [80]:

    
Ez = lambda x: 1/sqrt(0.3*(1+x)**3+0.7)

np.vectorize(Ez)
#Calculate comoving distance of a data point using the Redshift - This definition is based on the cosmology model we take. Here the distance for E-dS universe is considered. Also note that c/H0 ratio is cancelled in the equations and hence not taken.

def DC_LCDM(z):
  return integrate.quad(Ez, 0, z)



In [81]:

    
z=0.57
DMLCDM=DC_LCDM(z)
DMLC=z*(z/2+1)/(z+1)
print DMLCDM[0]
print DMLC









    



0.493380449608
0.46652866242



In [82]:

    
print DMLCDM[0]*4285.714/(1+z)
print DMLC*4285.714/(1+z)









    



1346.80732497
1273.50854773



In [83]:

    
th,phi=hu.randsphere(1000,'ang',theta_range=[0,0.25], phi_range=[0,0.25] )
thw,phiw=th,phi
#print th,phi
#print len(thw)



In [84]:

    
data = open("DAvspsi_z_0.57.dat",'w')
data.write("psi \t lfproper \t loproper \t DAf \t DAo \n")
#for i in range(len(thw)):
for thi,phii in zip(thw,phiw):
    for i in range(len(thw)-1):
    #print "i j"
        #print thw[i],phiw[i]
        #print thi,phii 
        #print cos(thw[i])*cos(thi)+sin(thw[i])*sin(thi)*cos(phiw[i]-phii)
        psi=acos(round(cos(thw[i])*cos(thi)+sin(thw[i])*sin(thi)*cos(phiw[i]-phii),6))
        xf1=DMLCDM[0]*sin(thw[i])*cos(phiw[i])
        yf1=DMLCDM[0]*sin(thw[i])*sin(phiw[i])
        zf1=DMLCDM[0]*cos(thw[i])
        xf2=DMLCDM[0]*sin(thi)*cos(phii)
        yf2=DMLCDM[0]*sin(thi)*sin(phii)
        zf2=DMLCDM[0]*cos(thi)
        lf=sqrt((xf1-xf2)**2+(yf1-yf2)**2+(zf1-zf2)**2)
        lfproper=lf/(1+z)
        xo1=DMLC*sin(thw[i])*cos(phiw[i])
        yo1=DMLC*sin(thw[i])*sin(phiw[i])
        zo1=DMLC*cos(thw[i])
        xo2=DMLC*sin(thi)*cos(phii)
        yo2=DMLC*sin(thi)*sin(phii)
        zo2=DMLC*cos(thi)
        lo=sqrt((xo1-xo2)**2+(yo1-yo2)**2+(zo1-zo2)**2)  
        loproper=lo/(1+z)
        if psi!=0:    
            DAf=lfproper/psi
            DAo=loproper/psi
            #print psi,DAf,DAo
            data.write("%f \t %f \t %f \t %f \t %f \n"%(psi,lfproper,loproper,DAf,DAo))
    thw=np.delete(thw,0)
    phiw=np.delete(phiw,0)
        #print thw
        #print phiw
    #print "array length"
    #print len(thw)
data.close()



In [85]:

    
DAvspsi=ascii.read("DAvspsi_z_0.57.dat")



In [86]:

    
print DAvspsi









    



  psi    lfproper loproper   DAf      DAo   
-------- -------- -------- -------- --------
0.201572 0.063238 0.059796 0.313724  0.29665
0.049523 0.015558 0.014711 0.314162 0.297064
0.106972 0.033601 0.031773 0.314115  0.29702
0.164783 0.051724 0.048909 0.313895 0.296811
0.187748 0.058914 0.055707 0.313791 0.296713
0.026982 0.008477 0.008015 0.314162 0.297064
0.027569 0.008666 0.008195 0.314346 0.297238
0.149652 0.046985 0.044428  0.31396 0.296873
0.021495 0.006759 0.006392 0.314474 0.297359
0.101414 0.031855 0.030121 0.314107 0.297012
     ...      ...      ...      ...      ...
0.032682  0.01027 0.009711 0.314253  0.29715
0.069785 0.021924 0.020731 0.314168  0.29707
0.040917 0.012859 0.012159  0.31426 0.297157
0.037605 0.011813  0.01117 0.314135 0.297039
0.047733 0.014996  0.01418 0.314166 0.297068
0.094162 0.029581 0.027971 0.314151 0.297054
0.085887 0.026983 0.025514 0.314169  0.29707
0.113516 0.035653 0.033713 0.314082 0.296988
0.008832 0.002774 0.002623 0.314041 0.296949
0.027786 0.008727 0.008252 0.314068 0.296975
0.035555 0.011171 0.010563 0.314186 0.297087
Length = 498297 rows



In [87]:

    
DAvspsi.sort('psi')



In [88]:

    
xint=DAvspsi['psi']
yint=1408*np.ones(len(xint))
#print xint



In [91]:

    
idx = np.argwhere(np.isclose(yint, DAvspsi['DAf']*4285.714 , atol=0.05)).reshape(-1)
print idx
print xint[idx]









    



[1699]
  psi   
--------
0.003162



In [92]:

    
print DAvspsi['DAo'][idx]*4285.714









    



     DAo     
-------------
1331.36991124



In [93]:

    
plt.plot(xint,yint,label="DA=1408")
plt.plot(xint[idx], yint[idx], 'ro')
#plt.plot(xint[idx]*np.ones(700),range(800,1500),'c',label="theta=inters")
plt.plot(xint[idx], DAvspsi['DAo'][idx]*4285.714, 'ro')
plt.plot(DAvspsi['psi'],DAvspsi['DAf']*4285.714,'g',label="DAf")
plt.plot(DAvspsi['psi'],DAvspsi['DAo']*4285.714,'r',label="DAo")
plt.legend(frameon=0, loc='upper right')
plt.show()



In [ ]: