User:Wesino/AgeOfUniverse

Age of Universe

Here is a code for computing the age of the universe times the Hubble time, as a function of the cosmological parameters Omega_matter and Omega_Lambda. It calculates the correction factor $F$ in the Age of the Universe article. This basically illustrates how increasing Lambda makes the universe older. In the code, we fix the radiation density, which is roughly equivalent to fixing the CMB temperature.

Following it is a Gnuplot script that will make a PostScript plot of this information, nicely labeled. On a Unix-type machine, if you have Gnuplot and Python_(language), then you can run:

$ python age.py > age.table
$ gnuplot age.gp

assuming you've called the Python program "age.py" and the Gnuplot script "age.gp" This creates a file called "age.eps" The Gnuplot script assumes that you've called the age data "age.table"

Python Program

Here's the Python program:

from math import *

def make_age_table():

    Omega_m_min = 0.01
    Omega_m_max = 1.2
    Omega_L_min = -0.2
    Omega_L_max = 1.0
    Nsteps = 50

    for omn in range(Nsteps):
        for oLn in range(Nsteps):
            Omega_m = Omega_m_min + omn*(Omega_m_max - Omega_m_min)/(Nsteps-1)
            Omega_L = Omega_L_min + oLn*(Omega_L_max - Omega_L_min)/(Nsteps-1)

            # WMAP matter-radiation equality redshift = 3454.  Use
            # this to fix the radiation Omega, using the fact that
            # WMAP also says that Omega_m = 0.266

            Omega_r0 = 0.266 / 3454.
            Omega_k = 1-Omega_r0-Omega_m-Omega_L
            
            age = calc_dimensionless_age( Omega_r0, Omega_m, Omega_k, Omega_L )

            print Omega_m, Omega_L, age
        print

def print_best_guess():
    """
    Print the best guess given various parameters determined by WMAP
    """

    print "Dimensionless age function assuming WMAP best fit of"
    print "Omega_m   0.266"
    print "Omega_L   0.732"
    print
    da = calc_dimensionless_age( 0.266/3454., 0.266, .002-0.266/3454., 0.732 )
    print "Yields: ", da
    print
    print "Best fit Hubble is 70 km/s/Mpc --> 13.97 Gyr"
    print "Combining yields age of ", da * 13.97, " Gyr"

def calc_dimensionless_age( Omega_r0, Omega_m0, Omega_k0, Omega_L0 ):
    """
    Compute the age function, so that the age of the universe is
    given by (1/H0)*F.  This computes F.
    """

    def f(a):
        return age_integrand( a, Omega_r0, Omega_m0, Omega_k0, Omega_L0 )

    return int_simp( f, 0., 1., 1001 )


def age_integrand( a, Omega_r0, Omega_m0, Omega_k0, Omega_L0 ):
    """
    This is the integrand in the age formula.

    It just scales the current values back to whenever a was.  This
    assumes that a=1 today.
    """

    I_denom = sqrt( Omega_r0 + a*(Omega_m0 + a*(Omega_k0 + a*a*Omega_L0)))

    return a/I_denom


def int_simp( f, a=0., b=1., nevals=10 ):
    """
    A simple Simpson's rule integrator.
    """

    # We need an odd number of evaluations.
    if nevals % 2 == 0:
        nevals += 1

    I = 0.
    dx = (b-a)/(nevals-1.)
    num_4s = (nevals-1)/2
    num_2s = (nevals-3)/2

    I = f(b) + f(a)

    # Add up the bits corresponding to the 4s.
    for j in range(num_4s):
        xj = a + dx*(2*j+1)
        I += 4*f(xj)

    # Now the bits corresponding to the 2s.
    for j in range(num_2s):
        xj = a + dx*2*(j+1)
        I += 2*f(xj)

    # Put it all together
    I *= dx/3.

    return I
   


if __name__ == "__main__":
    make_age_table()

Gnuplot Script

Here's the Gnuplot program

set terminal postscript
set output 'age.eps' 

set size square
set data style lines
set contour
set nosurface
set view 0,0
unset key
set grid
set xrange [0:1.2]
set yrange [-.2:1.0]
set cntrparam levels discrete  .666, .7, .8, .9, 1, 1.2, 1.5 

set label 1 "Age times current Hubble parameter" at .075,1.08

#set label 06 "0.6" at 1.35,-0.2
set label 0666 "0.667" rotate by 56 at 1.06,0.14
set label 07 "0.7" rotate by 60 at 1.07,0.50
set label 08 "0.8" rotate by 60 at 0.7,.73
set label 09 "0.9" rotate by 60 at 0.45,.82
set label 10 "1.0" rotate by 60 at 0.28, 0.85
set label 12 "1.2" rotate by 63 at 0.17,.87
set label 15 "1.5" rotate by 70 at 0.037, 0.80 

set label 100 "flat" rotate by -45 at .68,.38
set label 101 "closed" at .75,.50
set label 102 "open" at .55,.15
set label 103 "no CC" at .45,-0.04

set label 1000 "W" font "Symbol,24" at .5,-.40
set label 1001 "m" font "Times,18" at .60,-.43 

set label 1002 "W" font "Symbol,24" at 1.4,.4
set label 1003 "L" font "Symbol,18" at 1.50,.38

set label 200 at 0.266,0.732 point pt 4 ps 1
set label 201 at 1.,0 point pt 3 ps 1


splot 'age.table', 0.001*(x+y-1) +1, 0.001*y + 1

Enjoy!