In [363]:
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 [364]:
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 [365]:
z=0.35
DMLCDM=DC_LCDM(z)
DMLC=z*(z/2+1)/(z+1)
print DMLCDM[0]
print DMLC
In [366]:
print DMLCDM[0]*4285.714/(1+z)
print DMLC*4285.714/(1+z)
In [367]:
th,phi=hu.randsphere(1000,'ang',theta_range=[0,0.0025], phi_range=[0,0.0025] )
thw,phiw=th,phi
#print th,phi
#print len(thw)
In [368]:
data = open("DAvspsi_z_0.35.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),10))
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 [369]:
DAvspsi=ascii.read("DAvspsi_z_0.35.dat")
In [370]:
print DAvspsi
In [371]:
DAvspsi.sort('psi')
In [372]:
xint=DAvspsi['psi']*180/pi
yint=1050*np.ones(len(xint))
In [373]:
idx = np.argwhere(np.isclose(yint, DAvspsi['DAf']*4285.714 , atol=1)).reshape(-1)
print idx
print xint[idx]
In [374]:
print DAvspsi['DAo'][idx]*4285.714
In [376]:
plt.plot(xint,yint,label="DA=1050")
plt.plot(np.mean(xint[idx]), np.mean(yint[idx]), 'ro')
plt.plot(np.mean(xint[idx])*np.ones(500),range(800,1300),'c',label="theta=2.479")
plt.plot(np.mean(xint[idx]), np.mean(DAvspsi['DAo'][idx])*4285.714, 'ro')
plt.plot(DAvspsi['psi']*180/pi,DAvspsi['DAf']*4285.714,'g',label="DAf")
plt.plot(DAvspsi['psi']*180/pi,DAvspsi['DAo']*4285.714,'r',label="DAo")
plt.legend(frameon=0, loc='upper right')
plt.show()
In [375]: