#sample script to compute cosmic age for (flat or nonflat) w0wa model
# by Zhiqi Huang, huangzhq25@mail.sysu.edu.cn
#use 1/H0 as the time unit, but finally convert to Gyr for output
#if you want the result to be in h^{-1}Gyr, set H0 = 100 (i.e. h = 1)
import numpy as np
from scipy.integrate import quad
#density parameters
Omega_k = 0.0
Omega_r = 8.544e-5
Omega_m = 0.3
Omega_L = 1-Omega_k - Omega_r - Omega_m
w0 = -0.8
wa = 0.
z = 0. #where you want to compute the age of the universe
H0 = 70. #km/s/Mpc
## rho_DE/rho_{DE0} as a function of a
def rhoDE_ratio(a):
return np.exp(-3.*((1.+w0+wa)*np.log(a)+wa*(1.-a)))
## H as a function of a, in unit of H0
def Hofa(a):
return np.sqrt(Omega_r/a**4 + Omega_m/a**3 + Omega_k/a**2 + Omega_L*rhoDE_ratio(a))
## dt / da in unit of H0^{-1}
def dtbyda(a):
return 1.0/a/Hofa(a)
age, err = quad(dtbyda, 1.e-20, 1./(1.+z)) #age is in unit 1/H0, I start from 1.e-20 instead of 0 to avoid possible log(0) and 1/0 errors
#now convert to Gyr
Gyr = 1.e9*(365.2422*24*3600) #Gyr in SI unit
H0Gyr = H0*(1.e3/3.086e22)*Gyr #H0 in unit of Gyr^{-1}
print("At z = ", z, ", Age = " + str(age/H0Gyr) + " Gyr") #convert to Gyr