#sample script to compute cosmic age for LambdaCDM model 
# by Zhiqi Huang, for the course "General Relativity"
import numpy as np

#density parameters
Omega_k = 0.0
Omega_r = 0.0
Omega_m = 0.3
Omega_L = 1-Omega_k - Omega_r - Omega_m


##I use a_0=1 convention
Gyr = 1.e9*(365.2422*24*3600)  #Gyr in SI unit

H0Gyr = 70*(1.e3/3.086e22)*Gyr #H0 in unit of Gyr^{-1}


## H as a function of a, in unit of Gyr^{-1}
def Hofa(a):
    return H0Gyr*np.sqrt(Omega_r/a**4 + Omega_m/a**3 + Omega_k/a**2 + Omega_L)

## dt / da in unit of Gyr
def dtbyda(a):
    return 1.0/a/Hofa(a)


nsteps = 5000
astep = 1.0/nsteps
alist = np.arange(astep/2.0, 1.0, astep)
age = sum( [ dtbyda(a) for a in alist ] )*astep 
print("Age = " + str(age) + " Gyr")