Cosmology calculator

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 std model

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

Output section

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()