

In [21]:

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

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

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









    



1.31579183875
1.5202994012



In [24]:

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









    



1688.35554025
1950.76899039



In [25]:

    
th,phi=hu.randsphere(100,'ang')
thw,phiw=th,phi
#print th,phi
#print len(thw)



In [26]:

    
data = open("DAvspsi_z_2.34.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 [27]:

    
DAvspsi=ascii.read("DAvspsi_z_2.34.dat")



In [28]:

    
print DAvspsi









    



  psi    lfproper loproper   DAf      DAo   
-------- -------- -------- -------- --------
2.114055  0.68618  0.79283  0.32458 0.375028
0.701771 0.270824 0.312917 0.385916 0.445897
1.313443 0.481031 0.555795 0.366237 0.423159
0.381087 0.149223 0.172416 0.391571 0.452432
0.712472 0.274779 0.317487  0.38567 0.445613
0.743373 0.286155 0.330631 0.384941 0.444771
1.141039 0.425521 0.491657 0.372924 0.430885
 0.78772 0.302361 0.349355 0.383843 0.443502
1.024711 0.386253 0.446287 0.376938 0.435524
1.330487 0.486331  0.56192 0.365529 0.422341
     ...      ...      ...      ...      ...
2.459765 0.742555 0.857967 0.301881 0.348801
0.944206 0.358306 0.413995 0.379478 0.438459
 1.40022 0.507645 0.586546 0.362547 0.418896
 1.07348 0.402882   0.4655 0.375305 0.433637
1.412705 0.511397 0.590881 0.361998 0.418262
  0.8625 0.329347 0.380536 0.381852 0.441201
1.787403 0.614086  0.70953 0.343563 0.396961
0.535223 0.208343 0.240725 0.389265 0.449766
1.454242 0.523734 0.605135 0.360142 0.416117
0.517671 0.201667 0.233011 0.389565 0.450113
1.905759 0.642206 0.742022 0.336982 0.389358
Length = 4851 rows



In [29]:

    
DAvspsi.sort('psi')



In [30]:

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



In [37]:

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



In [38]:

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









    



     DAo     
-------------
1920.28272912



In [40]:

    
plt.plot(xint,yint,label="DA=1662")
plt.plot(np.mean(xint[idx]), yint[idx], 'ro')
plt.plot(xint[idx]*np.ones(1200),range(800,2000),'c',label="theta=inters")
plt.plot(xint[idx].mean(), 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 [33]: