

In [1]:

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

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

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









    



0.771427066428
0.75



In [19]:

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









    



1653.05788928
1607.14275



In [20]:

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



In [21]:

    
data = open("DAvspsi_z_1.0.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 [22]:

    
DAvspsi=ascii.read("DAvspsi_z_1.0.dat")



In [23]:

    
print DAvspsi









    



  psi    lfproper loproper   DAf      DAo   
-------- -------- -------- -------- --------
0.067866 0.026172 0.025445 0.385639 0.374928
0.010955 0.004212 0.004095 0.384542 0.373861
0.012806 0.004938   0.0048 0.385552 0.374843
0.043293 0.016694  0.01623 0.385595 0.374884
0.019596 0.007567 0.007357 0.386137 0.375412
0.024536 0.009467 0.009204 0.385839 0.375122
0.014283 0.005504 0.005351 0.385363 0.374659
0.007071 0.002713 0.002637 0.383617 0.372961
0.045895 0.017704 0.017212  0.38574 0.375026
0.040648 0.015676 0.015241 0.385664 0.374951
     ...      ...      ...      ...      ...
0.029632 0.011424 0.011107 0.385541 0.374832
0.081803 0.031546  0.03067 0.385632  0.37492
0.048025 0.018524 0.018009 0.385712 0.374998
   0.006 0.002319 0.002255 0.386551 0.375814
0.045567 0.017573 0.017085  0.38566 0.374947
0.033972 0.013098 0.012734 0.385551 0.374842
0.078875 0.030415  0.02957 0.385607 0.374896
0.036251  0.01398 0.013592 0.385657 0.374945
0.044926 0.017325 0.016843 0.385625 0.374914
0.003162 0.001275 0.001239  0.40307 0.391874
0.042735 0.016478  0.01602 0.385588 0.374878
Length = 495564 rows



In [24]:

    
DAvspsi.sort('psi')



In [25]:

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



In [36]:

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









    



[  6048  18573  27769 ..., 453061 453611 453723]
0.0402119643239



In [37]:

    
print np.mean(DAvspsi['DAo'][idx])*4285.714









    



1607.14233025



In [39]:

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