In [33]:
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 [52]:
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 [53]:
Ezo = lambda x: 1/sqrt(0.9*(1+x)**3+0.1*(1+x)**2)
np.vectorize(Ezo)
#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_OU(z):
return np.sinh(integrate.quad(Ezo, 0, z))
In [54]:
z=2.34
DMLCDM=DC_LCDM(z)
DMLC=DC_OU(z)
print DMLCDM[0]
print DMLC[0]
In [55]:
print DMLCDM[0]*4285.714/(1+z)
print DMLC[0]*4285.714/(1+z)
In [56]:
th,phi=hu.randsphere(1000,'ang',theta_range=[0,1], phi_range=[0,1] )
thw,phiw=th,phi
#print th,phi
#print len(thw)
In [57]:
data = open("DAvspsi_zou_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[0]*sin(thw[i])*cos(phiw[i])
yo1=DMLC[0]*sin(thw[i])*sin(phiw[i])
zo1=DMLC[0]*cos(thw[i])
xo2=DMLC[0]*sin(thi)*cos(phii)
yo2=DMLC[0]*sin(thi)*sin(phii)
zo2=DMLC[0]*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 [58]:
DAvspsi=ascii.read("DAvspsi_zou_2.34.dat")
In [59]:
print DAvspsi
In [60]:
DAvspsi.sort('psi')
In [61]:
xint=DAvspsi['psi']
yint=1662*np.ones(len(xint))
#print xint
In [62]:
idx = np.argwhere(np.isclose(yint, DAvspsi['DAf']*4285.714 , atol=0.05)).reshape(-1)
print idx
print xint[idx].mean()
In [63]:
print np.mean(DAvspsi['DAo'][idx])*4285.714
In [64]:
plt.plot(xint,yint,label="DA=1662")
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(xint[idx].mean(), 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 [64]:
In [64]: