Based on javascript code by Edward L. Wright
http://adsabs.harvard.edu/abs/2006PASP..118.1711W
Wright (2006, PASP, 118, 1711)
#!/usr/bin/env python
import sys
import numpy as np
from math import *
import matplotlib.pyplot as plt
from cProfile import label
class Universe:
def __init__(self): # Setup parameters and constants
self.H0 = 69.6 # Hubble constant
self.T0 = 2.72528 # CMB temperature in K
# initialize constants and parameters
self.WR = 0. # Omega(radiation)
self.WK = 0. # Omega curvaturve = 1-Omega(total)
self.WV = 0. # Omega vacuum
self.VWnu = 0. # Omega from massive neutrinos
self.mnue = 0.001 # electron neutrino mass ev
self.mnumu = 0.009 # mass of muon neutrino in eV
self.mnutau = 0.049 # mass of tau neutrino in eV
self.we = self.mnue/93 # Omega(nu(e))h^2
self.wmu = self.mnumu/93 # Omega(nu(mu))h^2
self.wtau = self.mnutau/93 # Omega(nu(tau))h^2
self.mnurel = 0.0005 # mass of neutrino that is just now relativistic in eV
self.c = 299792.458 # velocity of light in km/sec
self.Tyr = 977.8 # coefficent for converting 1/H into Gyr
self.DTT = 0.5 # time from z to now in units of 1/H0
self.DTT_Gyr = 0.0 # value of DTT in Gyr
self.age = 0.5 # age of Universe in units of 1/H0
self.age_Gyr = 0.0 # value of age in Gyr
self.zage = 0.1 # age of Universe at redshift z in units of 1/H0
self.zage_Gyr = 0.0 # value of zage in Gyr
self.DCMR = 0.0 # comoving radial distance in units of c/H0
self.DCMR_Mpc = 0.0
self.DCMR_Gyr = 0.0
self.DA = 0.0 # angular size distance
self.DA_Mpc = 0.0
self.DA_Gyr = 0.0
self.kpc_DA = 0.0
self.DL = 0.0 # luminosity distance
self.DL_Mpc = 0.0
self.DL_Gyr = 0.0 # DL in units of billions of light years
self.V_Gpc = 0.0
self.a = 1.0 # 1/(1+z), the scale factor of the Universe
self.az = 0.5 # 1/(1+z(object))
# Equation of state
# Note that the equation of state is given by
# w(a) = w + 2w'(1-a) following Linder (astro-ph/0402503).
self.w = -1 # equation of state, w = P/(rno*c^2)
self.wp = 0 # rate of change of equation of state, w(a) = w+2*wp*(1-a)
def nurho(self, mnu):
y = pow(1+pow(self. mnurel/mnu, 1.842), 1.0/1.842)
return y
def calcdcmt(self,WK,DCMR): # tangential comoving distance
ratio = 1.00
x = sqrt(abs(WK))*DCMR
if x > 0.1:
if WK > 0:
ratio = 0.5*(exp(x)-exp(-x))/x
else:
ratio = sin(x)/x
else:
y = x*x
if WK < 0: y = -y
ratio = 1. + y/6. + y*y/120.
return ratio*DCMR
def calcvcm(self,WK,DCMR): # comoving volume computation
ratio = 1.00
y=0
x = sqrt(abs(WK))*DCMR
if x > 0.1:
if WK > 0:
ratio = (0.125*(exp(2.*x)-exp(-2.*x))-x/2.)/(x*x*x/3.)
else:
ratio = (x/2. - sin(2.*x)/4.)/(x*x*x/3.)
else:
y = x*x
if WK < 0:
y = -y
ratio = 1. + y/5. + (2./105.)*y*y
return ratio*DCMR*DCMR*DCMR/3.
def correctannhilation(self,z): # correction for annihilations of particles not present now like e+/e-
lpz = log((1+1.0*z))/log(10.0)
dzage = 0
if lpz > 7.500 : dzage = 0.002 * (lpz - 7.500)
if lpz > 8.000 : dzage = 0.014 * (lpz - 8.000) + 0.001
if lpz > 8.500 : dzage = 0.040 * (lpz - 8.500) + 0.008
if lpz > 9.000 : dzage = 0.020 * (lpz - 9.000) + 0.028
if lpz > 9.500 : dzage = 0.019 * (lpz - 9.500) + 0.039
if lpz > 10.000 : dzage = 0.048
if lpz > 10.775 : dzage = 0.035 * (lpz - 10.775) + 0.048
if lpz > 11.851 : dzage = 0.069 * (lpz - 11.851) + 0.086
if lpz > 12.258 : dzage = 0.461 * (lpz - 12.258) + 0.114
if lpz > 12.382 : dzage = 0.024 * (lpz - 12.382) + 0.171
if lpz > 13.055 : dzage = 0.013 * (lpz - 13.055) + 0.188
if lpz > 14.081 : dzage = 0.013 * (lpz - 14.081) + 0.201
if lpz > 15.107 : dzage = 0.214
return dzage
def compute(self,z, H0=69.6, WM=0.286, WV=0.714): # computes the relevant values at the passed redshift
self.H0 = H0 # Hubble parameter
self.WM = WM # Omega(matter)
self.WV = WV # Omega(vacuum)
# Flat
# self.WV = 1.0-self.WM-0.2477 * pow(self.T0/2.72528,4)/(self.H0*self.H0) # Omega(vacuum) or lambda
# Open
#WV = 0.
# General
#WV = 0.714 # Entered value
self.h = self.H0/100.
self.WR = 2.477E-5 * pow(self.T0/2.7258,4)/(self.h*self.h) # does not include neutrinos, T0 = 2.72528
# rest mass omega*h^2 for the three neutrino types
self.we = (self.mnue/93.64) * pow(self.T0/2.72528,3)
self.wmu = (self.mnumu/93.90) * pow(self.T0/2.72528,3)
self.wtau = (self.mnutau/93.90) * pow(self.T0/2.72528,3)
# mass of nu that is just now relativistic
# evaluates at 3.151*kT with T = (4/11)^(1/3)*To and To=2.72528
# This is 6.13 K, and 1 eV is 11604.5 K
self.mnurel = 6.13*(self.T0/2.72528)/11604.5
self.Wnu = (self.we*self.nurho(self.mnue)+self.wmu*self.nurho(self.mnumu)+self.wtau*self.nurho(self.mnutau))/(self.h*self.h)
self.WK =1-self.WM-self.WR-self.WV
self.WM = self.WM-self.Wnu
self.az = 1.0/(1+1.0*z)
self.age = 0.
n=1000 # number of points in integrals
for i in range(n):
self.a = self.az*(i+0.5)/n
# rho(DE) = a^{-3-3*w_o-6*w'}*exp(6*w'*(a-1))*rho_o(DE)
# based on w = w_o+w_a*(1-a) with w_a = 2*w': Linder astro-ph/0402503
self.rhoV = self.WV*pow(self.a,-3-3*self.w-6*self.wp)*exp(6*self.wp*(self.a-1))
# get neutrino density corrected for kT/mc^2 by using lower mass
# instead of higher T:
self.Wnu = (self.we*self.nurho(self.mnue*self.a)+self.wmu*self.nurho(self.mnumu*self.a)+self.wtau*self.nurho(self.mnutau*self.a))/(self.h*self.h)
self.adot = sqrt(self.WK+((self.WM+self.Wnu)/self.a)+(self.WR/(self.a*self.a))+(self.rhoV*self.a*self.a))
self.age = self.age + 1./self.adot
self.zage = self.az*self.age/n
# correction for annihilations of particles not present now like e+/e-
# added 13-Aug-03 based on T_vs_t.f
dzage=self.correctannhilation(z)
self.zage = self.zage*pow(10.0,dzage)
self.zage_Gyr = (self.Tyr/self.H0)*self.zage
self.DTT = 0.0
self.DCMR = 0.0
# do integral over a=1/(1+z) from az to 1 in n steps, midpoint rule
for i in range(n):
self.a = self.az+(1-self.az)*(i+0.5)/n
self.rhoV = self.WV*pow(self.a,-3-3*self.w-6*self.wp)*exp(6*self.wp*(self.a-1))
self.Wnu = (self.we*self.nurho(self.mnue*self.a)+self.wmu*self.nurho(self.mnumu*self.a)+self.wtau*self.nurho(self.mnutau*self.a))/(self.h*self.h)
self.adot = sqrt(self.WK+((self.WM+self.Wnu)/self.a)+(self.WR/(self.a*self.a))+(self.rhoV*self.a*self.a))
self.DTT = self.DTT + 1./self.adot
self.DCMR = self.DCMR + 1./(self.a*self.adot)
self.DTT = (1.-self.az)*self.DTT/n
self.DCMR = (1.-self.az)*self.DCMR/n
self.age = self.DTT+self.zage
self.age_Gyr = self.age*(self.Tyr/self.H0)
self.DTT_Gyr = (self.Tyr/self.H0)*self.DTT
self.DCMR_Gyr = (self.Tyr/self.H0)*self.DCMR
self.DCMR_Mpc = (self.c/self.H0)*self.DCMR
# tangential comoving distance
self.DCMT = self.calcdcmt(self.WK, self.DCMR)
self.DA = self.az*self.DCMT
self.DA_Mpc = (self.c/self.H0)*self.DA
self.kpc_DA = self.DA_Mpc/206.264806
self.DA_Gyr = (self.Tyr/self.H0)*self.DA
self.DL = self.DA/(self.az*self.az)
self.DL_Mpc = (self.c/self.H0)*self.DL
self.DL_Gyr = (self.Tyr/self.H0)*self.DL
# comoving volume computation
self.VCM =self.calcvcm(self.WK, self.DCMR)
self.V_Gpc = 4.*pi*pow(0.001*self.c/self.H0,3)*self.VCM
# Output to the csv log file
# print('%1.1f,%1.3f,%1.3f' % (z,self.DCMR_Mpc, self.kpc_DA), file=f)
#return DCMR_Mpc
#return DA_Gyr, DTT_Gyr
return self.kpc_DA
def plot(Universe, xpoints, ypoints, axis, colour):
myLabel = r'$\Omega$m='+str(round(Universe.WM,3))
myLabel += ', $\Omega_\Lambda$='+str(round(Universe.WV,3))
myLabel += ', $\Omega$k='+str(round(Universe.WK,3))
myLabel += ', $\Omega$r='+str(round(Universe.WR,3))
myLabel += ', $\Omega$nu='+str(round(Universe.Wnu,3))
axis.plot(xpoints, ypoints, label=myLabel, linestyle='solid', color=colour)
axis.legend()
# Plot some charts
# Instantiate the object
Universe = Universe()
points = 50
i=-1
xpoints = np.zeros(points)
ypoints = np.zeros(points)
# Create plot and assign labels
fig, (ax1) = plt.subplots()
fig.suptitle(r'$\Lambda$CDM H$\_0$ = '+str(Universe.H0))
ax1.grid(True)
ax1.set_xlabel('redshift')
ax1.set_ylabel('Kpc/arcsec')
ax1.set_title('AngularScale(z)')
# Distance to farthest object https://en.wikipedia.org/wiki/GN-z11
zmax=15
step=zmax/points
for z in np.arange(0, zmax, step):
i=i+1
# get values for z =0 to 15
xpoints[i]=z
#compute(self,z, H0=69.6, WM=0.286, WV=0.714)
ypoints[i] = Universe.compute(z)
plot(Universe,xpoints,ypoints, ax1,'blue')
#ax2 = plt.twinx(ax1)
#ax2 = plt.twiny(ax1)
i=-1
xpoints = np.zeros(points)
ypoints = np.zeros(points)
for z in np.arange(0, zmax, step):
i=i+1
# get values for z =0 to 15
xpoints[i]=z
#compute(self,z, H0=69.6, WM=0.286, WV=0.714)
ypoints[i] = Universe.compute(z, 69.6, .1)
plot(Universe,xpoints,ypoints, ax1, 'red')
#ax3 = plt.twinx(ax2)
#ax3 = plt.twiny(ax2)
i=-1
xpoints = np.zeros(points)
ypoints = np.zeros(points)
for z in np.arange(0, zmax, step):
i=i+1
# get values for z =0 to 15
xpoints[i]=z
#compute(self,z, H0=69.6, WM=0.286, WV=0.714)
ypoints[i] = Universe.compute(z, 69.6, .1,.5)
plot(Universe,xpoints,ypoints, ax1, 'green')
plt.rcParams["figure.dpi"] = 120
plt.show()