

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

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









    



1.01893792017
1.05



In [5]:

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









    



1746.75060384
1799.99988



In [6]:

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

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

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



In [9]:

    
print DAvspsi









    



  psi    lfproper loproper   DAf      DAo   
-------- -------- -------- -------- --------
0.073732 0.030044 0.030959 0.407468 0.419889
0.099419 0.040503 0.041738   0.4074  0.41982
0.026534 0.010817 0.011146 0.407659 0.420086
0.101216 0.041235 0.042492 0.407397 0.419816
0.028179 0.011491 0.011841 0.407785 0.420216
0.115849 0.047192 0.048631 0.407361 0.419779
0.048626 0.019813 0.020417 0.407458 0.419879
0.049219 0.020062 0.020673   0.4076 0.420026
0.080508 0.032805 0.033805 0.407475 0.419896
0.055775 0.022726 0.023419 0.407464 0.419886
     ...      ...      ...      ...      ...
0.104919 0.042742 0.044045 0.407381   0.4198
0.140864 0.057365 0.059114 0.407237 0.419652
0.089091 0.036301 0.037408 0.407459 0.419881
0.102338 0.041692 0.042963 0.407392 0.419811
0.006481 0.002639  0.00272 0.407215 0.419629
0.053932 0.021978 0.022648 0.407504 0.419927
0.043247 0.017622 0.018159 0.407479   0.4199
0.137593 0.056036 0.057744 0.407258 0.419673
0.013491 0.005512  0.00568 0.408546    0.421
0.086699 0.035323   0.0364 0.407424 0.419844
0.100082 0.040775 0.042018 0.407415 0.419835
Length = 498300 rows



In [10]:

    
DAvspsi.sort('psi')



In [11]:

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



In [12]:

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









    



[  8993   9265   9714 ..., 247887 248240 248288]
0.0436125044771



In [13]:

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









    



1800.00111966



In [14]:

    
plt.plot(xint,yint,label="DA=1746.75060384")
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 [14]: