#sample script to compute cosmic age for (flat or nonflat) LambdaCDM 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.5e-5
Omega_m = 0.3
Omega_L = 1-Omega_k - Omega_r - Omega_m

z = 0. #where you want to compute the age of the universe
H0 = 70. #km/s/Mpc


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

## 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 divide-by-zero error
#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