The Museum of HP Calculators

General Relativity and Cosmology for the HP-41

Updated 12/12/05. New functions and Improvements noted in red.

This program is Copyright © 2005 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

 1°)  The Classical Tests

   a)   Advance of Periastron
   b)   Light Deflection
   c)   Gravitational Red-shift

 2°)  Gravitational Radiation on circular orbits

 3°)  Cosmology

   a)  Friedmann Universes:  Cosmological Constant = 0

     a-1)  Einstein-de Sitter model  q0 = 1/2     ( q = deceleration parameter )
     a-2)  Spherical models  q0 > 1/2
     a-3)  Hyperbolic models  0 < q0 < 1/2

   b)  Empty Universes:   L0 + q0 = 0

     b-1)  Milne model  q0 = 0
     b-2)  de Sitter model  q0 = -1
     b-3)  -1 < q0 < 0
     b-4)   q0 > 0
     b-5)   q0 < -1

   c) Eddington-Lemaitre Universes ( critical parameters )

   d) Euclidean Universes

     d-1)   -1 < q0 < 1/2
     d-2)    q0 > 1/2

   e)  More general models

     e-1) Numerical Integration ( program#1 )
     e-2) Numerical Integration ( program#2 )
     e-3) Elliptic Integrals ( Big Bang Universes only )

   f)  Non-vanishing Pressure   ( new )

     f-1) Pure Radiation  ( new )
     f-2) Radiation+Dust:  Numerical Integration( new )
     f-3) Radiation+Dust:  Elliptic Integrals( new )
 

Astronomical Data:

              Constant of gravitation = G = 6.67259 10-11  m3 kg-1 s-2
              Heliocentric Gravitational Constant = 1.3271244 1020  m3 s-2
              Geocentric Gravitational Constant = 3.986004356 1014  m3 s-2
              Mass of the Sun = M = 1.98892 1030  kg
              Mass of the Earth = m = 5.97370 1024  kg
              Solar Radius = 696265 km
              Earth Equatorial Radius = 6378137 m
              Speed of light = c = 299792458 m/s
              Astronomical Unit = 1AU = 149597870.6 km
              1 Light-year = 9.4607304726  1015 m
              1 Megaparsec = 1 Mpc = 3.085677581 1022 m
              Stefan's Constant = a = 7.565656  10 -16 J.m3.K -1
              Temperature of the Cosmic Background = T = 2.736 K
              Hubble Constant = H0 ~ 71 km/s/Mpc ~ 2.301 10-18 s

-In the following programs, I've choosen  1/H0 = 13.7 109 years which corresponds to  H0 = 71.37169499 km/s/Mpc = 2.312999111 10-18 second
-These decimals are only given to facilitate comparisons: Hubble Constant is known very inaccurately!
-The mean density of matter is of the order of  ( rho )mat ~  3  10 -28 kg/m3 , perhaps 10 or 100 times more!
  ( possible dark matter , non-baryonic or exotic matter ... )

-Some of these programs use synthetic registers M , N , O , P , Q   which may of course be replaced by any unused data registers.
-Integrals are denoted  §
 

1°) The Classical Tests
 

        a) Relativistic Advance of Periastron
 

-We assume a celestial body of mass M creates a spherically symmetric gravitational field.
-In this case, the metric has the following expression:

   ds2 = (1-2GM/(c2.r)) c2.dt2 - dr2/(1-2GM/(c2.r)) - r2(d(theta)2 + sin2(theta) d(phi)2 )     Schwarzschild metric

 ( r , theta , phi  = spherical coordinates )

-Provided the gravitational field is not too strong, the orbit of a planet is a slowly rotating ellipse.
-The 3 following short routines use the formulae:

        Advance of periastron = 6.PI.GM/(a.c2(1-e2)) radians per revolution
        Third Kepler's law:  a3/T2 = GM/(4.PI2)

 where    M = Mass of the system  ,  T = orbital period  ,  a = semimajor axis  ,  e = eccentricity.
 

01  LBL "MTE"
02  X^2
03  SIGN
04  LASTX
05  -
06  .198956
07  X<>Y
08  /
09  RDN
10  ST/ T
11  /
12  X^2
13  3
14  1/X
15  Y^X
16  R^
17  *
18  END

( 33 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             M             /
           Y             T            M
           X             e         adv(p)

  where    M = Mass ( expressed in Solar Masses )
                T = Orbital Revolution ( in days )
                e = eccentricity
and  adv(p) = the advance of periastron in degrees per year
 

Example1:    The binary Pulsar  PSR J0737-3039  has the following characteristics:

    m1 = 1.25 Solar Mass  ;   m2 = 1.34 Solar Mass   whence   M = m1 +  m2 = 2.59
    T  = 2.45 hours
    e  =  0.09

    2.59  ENTER^
    2.45  ENTER^  24  /
    0.09  XEQ "MTE"  >>>>  adv(p) = 16.97 deg/year

Example2:    For Mercury,  M = Mass of the Sun = 1
                                             T = 87.969 days
                                              e = 0.2056

        1       ENTER^
    87.969  ENTER^
    0.2056     R/S       >>>>   0.00011939 deg/year.      Multiplying this value by  36 E4  it yields:   42.98 arcseconds/century

-The advance of perihelion for Mercury is the first classical test of general relativity.
-Of course, Newtonian perturbations caused by other planets are to be added to this relativistic effect.
-This routine assumes the correction is small.
-If  M is expressed in kg and  T in seconds, replace line 06 with  2.124235 E-13
-Similar programs may be written if we know  a , T , e  instead of M , T , e:
 

01  LBL "ATE"
02  X^2
03  SIGN
04  LASTX
05  -
06  519.464
07  X<>Y
08  /
09  ST* Z
10  CLX
11  3
12  Y^X
13  /
14  *
15  END

( 30 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             a            a
           Y             T            a
           X             e         adv(p)

where   a = Semimajor axis ( in Astronomical Units )             If  a is expressed in meters and T in seconds,
           T = Orbital Revolution ( in days )                              replace line 06 by  1.497084 E-5
            e = eccentricity
 

-And if we know  M , a , e:
 

01  LBL "MAE"
02  X^2
03  SIGN
04  LASTX
05  -
06  93808
07  *
08  RCL Y
09  X^2
10  *
11  ST/ T
12  RDN
13  /
14  SQRT
15  *
16  END

( 30 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             M            /
           Y              a            /
           X              e         adv(p)

where  M = Mass ( expressed in Solar Masses )                   If  M is expressed in kilograms and a in meters
            T = Orbital Revolution ( in days )                              replace line 06 with  3.039848 E22
            e = eccentricity
 

        b) Light Deflection

-If a photon passes at a distance d of a celestial body of mass M,
  its trajectory is deflected by an angle  LD = 4.GM/(c2d)  radians.
-The (very ) short following routine computes LD in seconds of arc.
 

01  LBL "LD"
02  /
03  1.7498
04  *
05  END

( 17 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Y              M            /
           X              d         LD(")

 where  M is expressed in Solar Masses
   and    d  --------------------  radii.

Example:   with  M = 1 for the Sun and d = 1 it yields of course  1.75"

-This effect of light deflection was first observed at the Solar eclipse of 1919.
-If M is expressed in kg and d is expressed in meters, replace line 03  by  6.1255 E-22
-This program is only valid for small deflections.
-Assuming spherical symmetry, the exact differential equation is:

  d2u/d(phi)2 + u = 3.GM/c2 u2     where   u = 1/r    and ( r, phi ) = polar coordinates.

-A "funny" solution is  r = Constant = 3.GM/c2 = 1.5 * Schwarzschild radius
-Thus, a photon may have a circular trajectory !
 

        c) Gravitational Red-Shift
 

-We suppose the source and the observer are static in a spherically symmetric gravitational field created by a celestial body of mass M.
-In this case,  z = d(lambda)/lambda = ((1-2GM/c2/rO)/(1-2GM/c2/rS))1/2  - 1      ( r = distance from the origin )
 

01  LBL "GRZ"
02  4.2416 E-6
03  R^
04  *
05  CHS
06  ENTER^
07  R^
08  /
09  X<> Z
10  /
11  RCL X
12  RCL Z
13  -
14  X<> Z
15  1
16  ST+ Z
17  +
18  ST* Y
19  X<>Y
20  SQRT
21  +
22  /
23  END

( 45 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             M            a
           Y             rS            a
           X             rO         adv(p)

where    M  = Mass of the gravitational body           ( in Solar mass )
              rS = distance source-gravitational body     ( in Solar radius )
              rO = distance observer-gravitational body  ----------------

Example:     M = 2 ,  rS = 3 ,  rO = 4

    2  ENTER^
    3  ENTER^
    4  XEQ "GRZ"   >>>>  z = 3.5347 10-7

-If  M is expressed in kg and  rS ,  rO are in meters, replace line 02 by  14849 E-31
-When M = mass of the Earth, rO = radius of the Earth and rO-rS = 22.496 m, it yields  z ~  2.46 10-15
-This very small effect has been measured and confirmed by Pound & Rebka in 1960

-If both GM/(c2r) are small (  GM/(c2r) << 1 )  we can use a much shorter routine:

01  LBL "GRZ2"
02  1/X
03  CHS
04  X<>Y
05  1/X
06  +
07  *
08  2.1208 E-6        or  74243 E-32   if   M is expressed in kg  and  r in meter
09  *
10  END
 

2°) Gravitational Radiation on Circular Orbits
 

-We now assume that 2 celestial bodies move on cicular orbits around their common center of gravity.
-Because of gravitational radiation, the system looses a small energy, and their trajectories are actually spirals.
-The distance r between these 2 masses ( m1 & m2 ) is slowly decreasing according to the differential equation:

   dr/dt = -64G3 m1.m2 ( m1 + m2 ) / ( 5.c5.r3 )

-If  r = r0  at an instant  t = 0  the following program calculates at what time T we'll have  r = r1    (  r1 <  r0 )

01  LBL "SPI"
02  X^2
03  X^2
04  X<>Y
05  X^2
06  X^2
07  -
08  ABS
09  RDN
10  ST* Z
11  +
12  *
13  /
14  32114 E13
15  *
16  END

( 32 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T             m1            m1
           Z             m2            m1
           Y              r0            m1
           X              r1            T

  where  mi are expressed in Solar Masses and  rj  are expressed in Astronomical Units
             The result T is expressed in years

Example:   Let's take once again the binary pulsar PSR J0737-3039      m1 = 1.25  ,   m2 = 1.34    ( we neglect the small eccentricity )
 -Currently,  r is approximately  0.00586 AU
 -Find in how many time we'll have r = 0.001 AU

      1.25     ENTER^
      1.34     ENTER^
   0.00586  ENTER^
    0.001     XEQ "SPI"  >>>>  87 Million years

-The coalescence time is estimated  85 Million years ( it's the shortest currently known ).
-If mi are in kg and  rj are in meters, replace line 14 with  1.592026 E71
 

3°) Cosmology
 

-Einstein's equations   Sij - (1/2) gij ( S + 2.Lambda ) = ( 8.PI.G/c4 ) Tij   are now applied to a homogeneous and isotropic Universe.
    (  Sij is the tensor of curvature , S = SijSij , Lambda is the cosmological constant and  Tij describes the content of the Universe  )

-This "cosmological principle" leads to the Robertson-Walker metric:  ds2 = c2dt2 - R2(t)/(1-kr2/4)2 . ( dr2 + r2d(theta)2 + r2sin2(theta)d(phi)2 )

      r , theta , phi  are the comoving spherical coordinates
    R(t) is the scale factor = the "radius" of the Universe
     k = +1  for Spherical   Universes
     k =  0   for Euclidean   ---------
     k = -1  for Hyperbolic ---------

-Usually,  Tij = ( (rho).c2 + p ) uiuj - p gij      ( gij is the metric tensor ; rho = density of matter and energy ; p = pressure ; ui = dxi/ds = 4-velocity )
  and Einstein's equations yield:

  2 R(t).d2R/dt2 + (dR/dt)2 + k.c2 = ( -(8.PI.G/c2 ).p + (Lambda).c2 ).R2(t)
     (dR/dt)2 + k.c2  = ( (8.PI.G/3) (rho) + (Lambda/3).c2 ).R2(t)

_________________________________________________________________________________________________________________________

-If R(t) is constant and p = 0, this system leads to Einstein's static model:    k = +1 , R = RE = c.(4.PI.G.(rho))-1/2 ; Lambda = 4.PI.G.(rho)/c2 = 1/R2

   With  rho = 5 10 -28 kg/m3 it yields   R = 4.63 1026 m ~ 15 Gpc  ,  Lambda = 4.665 10 -54 m-2
   Einstein's Universe is a 3-sphere and its volume is  2.PI2.R3 = 1.959 1081 m3 ~ 66700 Gpc3
   Its mass is M = 9.796 1053 kg

_________________________________________________________________________________________________________________________

-To simplify the equations when R is not constant, we now assume p = 0 ( vanishing pressure ) and we use the following parameters:

           H = (1/R).dR/dt                   =  Hubble constant
           q  = -R.(d2R/dt2)/(dR/dt)2   =  deceleration parameter
  (Omega)mat = 8.PI.G.(rho)/(3H2) =  density parameter
  (Omega)vac = (Lambda.c2)/(3H2) =  L

-The equations reduce to:

          2.q + 2.L = 8.PI.G.(rho)/(3H2)
     2.q + 3.L - 1 = k.c2/(H2R2)

  and lead to the integral      H0.t = §  [ ( 2q0+2.L0 )/y + ( 1-2q0-3L0 ) + L0.y2 ] -1/2 dy

   where   y = R/R0  , the subscript 0  meaning a current value

-Knowing the cosmological redshift z of a galaxy, and the values of the current parameters q0 , L0
 the following programs computes several "distances", the recessional velocity and the age of the Universe.

            D   =  light-travel time distance ( in light-years ) = light-travel time ( in years )
             D0 = comoving distance ( in light-years )
             DL = luminosity-distance ( in light-years )
             VR = recessional velocity ( speed of light = 1 )
              t0  = age of the Universe ( in years )

     D  = (c/H0)  §y(em)1  [ ( 2q0+2.L0 )/y + ( 1-2q0-3L0 ) + L0.y2 ] -1/2 dy               y(em) = y at the instant of emission     z + 1 = R0/R = 1/yem
     D0 = (c/H0)  §y(em)1  [ ( 2q0+2.L0 ).y + ( 1-2q0-3L0 ).y2 + L0.y4 ] -1/2 dy           D0 = Dnow

                                                                F(x) = Sinh(x)  if  k = -1
     DL = R0 ( z + 1 ) F(D0/R0)    where    F(x) =   x         if  k = 0
                                                               F(x) = Sin(x)    if  k = +1

     VR = H0.D0
      t0  =  (c/H0)  §01  [ ( 2q0+2.L0 )/y + ( 1-2q0-3L0 ) + L0.y2 ] -1/2 dy

-We also have  k = sign( 2q0+3L0-1 )   and   R0 = (c/H0) k.( 2q0+3L0-1 ) -1/2   R0 may be choosen arbitrarily ( positive ) if k = 0

-Since  z + 1 = R0/R = 1/y , these programs may also be used to draw the graph of the function R(t)
    z > 0  for the past ,  -1 < z < 0  for the future.

-Here is a summary of the different models:
 

_______________________________________________MANY UNIVERSES___________________________________________________
 

-According to the values of q and L, the different kinds of Universes may be represented as follows:
 

  V    Hyp    E       Sph        q                C1
     *                e                   |            c
         *               e     Dom2 |         c
             *             e              |      c
                 *  Dom1e           |    c
                     *           e        |  c
                         *          e   1| c
                             *         e  |c                          Dom4
                                 *        eA(0;1/2)                                                A = Einstein-de Sitter model
                                     *     | e
     ------------------------O|D3e--------------------------------------- L
                                            |  *   e       1              2
                                            |      *  e
                                            |          * e
                                       -1 |                cB(1;-1)                                 B = de Sitter model
                                            |                  * c
                                            |                      *   c
                                            |                          *       c
                                       -2 |                               *            c
                                            |                                  *   Dom5       c
                                            |                                       *                           C2
                                                                                        W

-Line (VW): q + L = 0         represents empty Universes. On the left of this line, the density of matter is negative:  q + L < 0
-Line (BE): 2q + 3L -1 = 0  ----------  Euclidean Universes.  On the left of this line, Universes are hyperbolic. On the right, they are spherical.
 
 

-Domains1&2:  VOAE  & EAC1  there is a Big Bang and a Big Crunch

       R(t)
           |                *
           |      *                 *
           |  *                          *
   --------|*------------------------------*---------  t
          O

-Domains3&4:  OAB  & C1ABC2  there is a Big Bang and no Big Crunch. R(t) has a point of inflexion and R(t) tends to infinity as t tends to infinty
                          Our Universe seems to be inside triangle OAB, near the midpoint of [OB]

       R(t)
           |                                     *
           |                                  *
           |                             *
           |                      *
           |           *
           |    *
     ------|*---------------------------------------------  t
          O

-Domain5:    WBC2  there is no Big Bang and no Big Crunch. R(t) has a positive minimum ( bounce models )

       R(t)
           |
           |  *                                *
           |   *                             *
           |     *                        *
           |          *              *
           |                  *
     ------|-----------------------------------  t
          O

-Critical parameters:

-Lines  )AcccC1) and )BcccC2)  represent "critical" cosmological constants. On theses lines, we have:  ( 3L+2q-1)3 = 27 L(q+L)2
-The straight-line  R = RE is an asymptote of the curve R = R(t) in both cases:

       R(t)
          |                                                                  *
           |                                                         *               Einstein-Lemaitre models:   Line )BC2)      There is no Big Bang and  R(t)  = RE  for t = - infinity
           |                                           *
           |                          *
           |*
           |---------------------------------------------------------------------------- Einstein's model   R(t) = RE = constant
           |                                                                       *
           |                           *
           |           *
           |    *                                                                    Einstein-Eddington models:  Line )AC1)      There is a Big Bang and R(t)  = RE  for t = + infinity
     ------|*-------------------------------------------------------------------------- t
 
 

_________________________________________________________________________________________________________________________
 
 

        a) Friedmann Universes: No Cosmological Constant , L0 = 0

         a-1)  Einstein-de Sitter Model:  q0 = 1/2 ; L0 = 0

-This is an Euclidean Universe:   2q0 + 3L0 - 1 = 0     whence  k = 0

Formulae:

    D = (c/H0).(2/3).( 1 - (1 + z )-3/2 )
   D0 = 2(c/H0).( 1 - (1 + z )-1/2 )
   DL = D0 ( 1 + z )                               R(t) = (3H0t/2)2/3    if we choose  R0 = 1
   VR = H0 D0
    t0 =  (2/3)H0
 

01  LBL "EdS"
02  1
03  +
04  ENTER^
05  ENTER^
06  SQRT
07  1/X
08  CHS
09  ENTER^
10  R^
11  /
12  1
13  ST+ Z
14  +
15  STO M
16  RDN
17  STO Z
18  ST+ Z
19  ST* Y
20  137 E8
21  ST+ X
22  ST* M
23  ST* Z
24  *
25  3
26  ST/ M
27  ST/ L
28  CLX
29  X<> M
30  END

( 53 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y              /            D0
           X              z            D
           L              /            t0

 

Example:      z = 7

    7  XEQ "EdS"  >>>>  D  =   8.730 109 ly
                            RDN   D0 = 17.713 109 ly
                            RDN   DL = 141.701 109 ly
                            RDN   VR =  1.2929
                         LASTX  t0  =  9.133 109 years
 

         a-2)  Spherical Models:  q0 > 1/2 ; L0 = 0

-Here,  k = +1

Formulae:

    D = (c/H0/(2q0-1)).[ -1 + (2q0z+1)1/2/(z+1) - (q0/(2q0-1)1/2) ( Arccos ((2q0-1)/q0-1) - Arccos ((2q0-1)/(q0(z+1)) -1) ]
   D0 = (c/H0/(2q0-1)1/2).[ - Arccos ((2q0-1)/q0-1) + Arccos ((2q0-1)/(q0(z+1)) -1) ]
   DL = R0 ( 1 + z ) Sin D0/R0                     R0 = c/H0/(2q0-1)1/2
   VR = H0 D0
    t0 =  (1/H0/(2q0-1)) [ -1 - (q0/(2q0-1)1/2) ( Arccos ((2q0-1)/q0-1) - PI ]

-The period of these Universes ( time between a Big Bang and a Big Crunch ) is  T = (1/H0)(2.pi.q0).(2q0-1) -3/2
 

01  LBL "FRSPH"
02  DEG
03  1
04  +
05  STO M
06  X<>Y
07  ENTER^
08  STO P        ( synthetic )
09  ST+ X
10  LASTX
11  -
12  STO N
13  X<>Y
14  /
15  ENTER^
16  ENTER^
17  SIGN
18  -
19  ACOS
20  CHS
21  STO T
22  X<> Z
23  /
24  PI
25  R-D
26  ST+ T
27  SIGN
28  -
29  ACOS
30  +
31  STO O
32  RCL P
33  RCL N
34  SQRT
35  /
36  D-R
37  ST* Z
38  *
39  RCL M
40  RCL P
41  ST+ X
42  ST* Y
43  -
44  1
45  ST- T
46  +
47  SQRT
48  RCL M
49  ST- Y
50  /
51  +
52  X<> M
53  RCL O
54  SIN
55  *
56  RCL N
57  ST/ M
58  ST/ Z
59  SQRT
60  ST/ Y
61  R-D
62  ST/ O
63  X<> Z
64  SIGN
65  X<> O
66  STO Z
67  137 E8
68  ST* L
69  ST* M
70  ST* Y
71  ST* Z
72  X<> M
73  CLA
74  END

( 121 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z             /            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:      q0 = 2  ;  z = 3

    2  ENTER^
    3  XEQ "FRSPH"  >>>>  D  =  5.871 109 ly
                                 RDN   D0 =  9.482 109 ly
                                 RDN   DL = 29.474 109 ly
                                 RDN   VR =  0.6921
                               LASTX  t0  =  6.477 109 years

Note:  For  q0 = 0.50001  the results are still approximately correct ( almost those of Einstein-de Sitter model )
      but with q0 = 0.500001 ,  D and t0 are wrong! ( due to roundoff errors )
 

         a-3)  Hyperbolic Models:  0 < q0 < 1/2 ; L0 = 0

-Here,  k = -1

Formulae:

    D = (c/H0/(1-2q0)).[ 1 - (2q0z+1)1/2/(z+1) - (q0/(1-2q0)1/2) ( Arccosh ((1-q0)/q0) - Arccos ((1-2q0)/(q0(z+1)) +1) ]
   D0 = (c/H0/(1-2q0)1/2).[  Arccosh ((1-q0)/q0) - Arccosh ((1-2q0)/(q0(z+1)) +1) ]
   DL = R0 ( 1 + z ) Sinh D0/R0                               R0 = c/H0/(1-2q0)1/2
   VR = H0 D0
    t0 =  (1/H0/(1-2q0)) [ 1 - (q0/(1-2q0)1/2)  Arccosh ((1-q0)/q0)  ]
 

01  LBL "FRHYP"
02  1
03  +
04  STO M
05  X<>Y
06  ENTER^
07  STO P         ( synthetic )
08  ST+ X
09  SIGN
10  LASTX
11  -
12  STO N
13  X<>Y
14  /
15  ENTER^
16  ENTER^
17  X^2
18  LASTX
19  ST+ X
20  +
21  SQRT
22  +
23  LN1+X
24  X<>Y
25  RCL M
26  /
27  ENTER^
28  X^2
29  LASTX
30  ST+ X
31  +
32  SQRT
33  +
34  LN1+X
35  -
36  STO O
37  RCL P
38  RCL N
39  SQRT
40  /
41  CHS
42  ST* Z
43  *
44  RCL M
45  RCL P
46  ST+ X
47  ST* Y
48  -
49  1
50  ST+ T
51  +
52  SQRT
53  RCL M
54  ST- Y
55  /
56  -
57  X<> M
58  RCL O
59  E^X
60  ENTER^
61  1/X
62  -
63  *
64  2
65  /
66  RCL N
67  ST/ M
68  ST/ Z
69  SQRT
70  ST/ Y
71  ST/ O
72  X<> Z
73  SIGN
74  X<> O
75  STO Z
76  137 E8
77  ST* L
78  ST* M
79  ST* Y
80  ST* Z
81  X<> M
82  CLA
83  END

( 130 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:      q0 = 0.4  ;  z = 7

   0.4  ENTER^
    7  XEQ "FRHYP"  >>>>  D  =   9.087 109 ly
                                 RDN   D0 = 18.708 109 ly
                                 RDN   DL = 159.140 109 ly
                                 RDN   VR =  1.3655
                               LASTX  t0  =  9.534 109 years
 

        b) Empty Universes:   q0 + L0 = 0

         b-1)  Milne Model:   q0 = 0

  k = -1  ( hyperbolic geometry )

Formulae:

    D = (c/H0).z/(z+1)
   D0 = 2(c/H0).Ln(z+1)
   DL = R0 z(z+2)/2           R0 = c/H0        R(t) = c.t
   VR = H0 D0
    t0 =  1/H0
 

01  LBL "MLN"
02  ENTER^
03  STO M
04  1
05  +
06  ST/ M
07  LN
08  X<>Y
09  1
10  LASTX
11  +
12  *
13  2
14  /
15  X<>Y
16  137 E8
17  ST* M
18  ST* Z
19  *
20  0
21  X<> M
22  END

( 39 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y              /            D0
           X              z            D
           L              /            t0

Example:        z = 7

    7  XEQ "MLN"  >>>>  D  =  11.988 109 ly
                              RDN   D0 = 28.488 109 ly
                              RDN   DL = 431.550 109 ly
                              RDN   VR =  2.0794
                             LASTX  t0  = 13.700 109 years

         b-2) de Sitter  Model:   q0 = -1

  k = 0  ( Euclidean geometry )

Formulae:

    D = (c/H0).Ln(z+1)
   D0 = (c/H0).z
   DL = (c/H0).z(z+1)
   VR = H0 D0
    t0 =  infinity              R(t) = exp(H0.t)      the mimimum = 0  when  t = minus infinity
 

01  LBL "dST"
02  ENTER^
03  STO Z
04  1
05  +
06  ST* Z
07  LN
08  STO M
09  CLX
10  137 E8
11  ST* M
12  ST* Z
13  *
14  DEG
15  90
16  TAN
17  SIGN
18  CLX
19  X<> M
20  END

( 39 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y              /            D0
           X              z            D
           L              /            t0

Example:        z = 7

    7  XEQ "dST"  >>>>  D  =  28.488 109 ly
                           RDN   D0 =  95.900 109 ly
                           RDN   DL = 767.200 109 ly
                           RDN   VR =  7.0000
                         LASTX  t0  =  9.9999 1099 years = infinity

         b-3)  -1  <  q0  <  0

-Recent observations suggest that  q0 is of the order of  -1/2

Formulae:       k = -1   ( hyperbolic geometry )

    D = (c/H0/(-q0)1/2). [ Arcsinh (1/b) - Arcsinh (1/b/(z+1)) ]      where   b = ((1+q0)/(-q0))1/2
   D0 = (c/H0/(1+q0)1/2). Ln ( ( 1/b2 + (z+1)2 )1/2 + z+1 )/((1/b2+1)1/2+1 )
   DL = R0 ( 1 + z ) Sinh D0/R0          R0 = c/H0/(1+q0)1/2       R(t) = R0.b Sinh ( H0.t.(-q0)1/2 )
   VR = H0 D0
    t0 =  1/H0/(-q0)1/2 Arcsinh 1/b
 

01  LBL "VAC"
02  1
03  +
04  STO M
05  X<>Y
06  ENTER^
07  CHS
08  STO N
09  LASTX
10  -
11  /
12  STO O
13  SQRT
14  LASTX
15  1
16  STO Q
17  +
18  SQRT
19  ST+ Q
20  +
21  ENTER^
22  R^
23  ST* Y
24  X^2
25  RCL O
26  +
27  SQRT
28  STO P         ( synthetic )
29  RCL O
30  SQRT
31  +
32  /
33  LN
34  X<> M
35  ENTER^
36  X<> P
37  +
38  RCL Q
39  /
40  LN
41  STO O
42  X<> L
43  ENTER^
44  1/X
45  -
46  RCL P
47  *
48  2
49  /
50  X<>Y
51  LN
52  RCL N
53  SQRT
54  ST/ M
55  /
56  1
57  RCL N
58  -
59  SQRT
60  ST/ O
61  ST/ Z
62  RDN
63  SIGN
64  X<> O
65  STO Z
66  137 E8
67  ST* L
68  ST* M
69  ST* Y
70  ST* Z
71  X<> M
72  CLA
73  END

( 112 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:      q0 = -0.3 ;  z = 7

 -0.3  ENTER^
    7   XEQ "VAC"   >>>>  D  =  13.341 109 ly
                                RDN   D0 =  32.552 109 ly
                                RDN   DL = 469.215 109 ly
                                RDN   VR =  2.3761
                              LASTX  t0  =  15.386 109 years

         b-4)  q0  >  0

   k = -1   ( hyperbolic geometry )

Formulae:

    D = (c/H0/q01/2). [ Arcsin (1/b) - Arcsin (1/b/(z+1)) ]      where   b = ((1+q0)/q0)1/2
   D0 = (c/H0/(1+q0)1/2). Ln ( ( ( -1/b2 + (z+1)2 )1/2 + z+1 )/((-1/b2+1)1/2+1 )
   DL = R0 ( 1 + z ) Sinh D0/R0          R0 = c/H0/(1+q0)1/2       R(t) = R0.b Sin ( H0.t.q01/2 )
   VR = H0 D0
    t0 =  1/H0/q01/2 Arcsin 1/b
 

01  LBL "VAC2"
02  DEG
03  1
04  +
05  X<>Y
06  ENTER^
07  STO N
08  LASTX
09  +
10  /
11  STO M
12  CHS
13  X<>Y
14  X^2
15  +
16  SQRT
17  +
18  1
19  RCL M
20  -
21  SQRT
22  1
23  +
24  /
25  LN
26  STO O
27  X<> L
28  ENTER^
29  1/X
30  -
31  *
32  2
33  /
34  RCL M
35  SQRT
36  ST/ Z
37  ASIN
38  ENTER^
39  R^
40  1/X
41  ASIN
42  -
43  D-R
44  STO M
45  RDN
46  D-R
47  RCL N
48  SQRT
49  ST/ M
50  /
51  1
52  RCL N
53  +
54  SQRT
55  ST/ O
56  ST/ Z
57  RDN
58  SIGN
59  X<> O
60  STO Z
61  137 E8
62  ST* L
63  ST* M
64  ST* Y
65  ST* Z
66  X<> M
67  CLA
68  END

( 101 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:      q0 =  0.3 ;  z = 7

   0.3  ENTER^
    7   XEQ "VAC2"   >>>>  D  =  11.031 109 ly
                                 RDN   D0 =  25.737 109 ly
                                 RDN   DL = 403.673 109 ly
                                 RDN   VR =  1.8786
                               LASTX  t0  =  12.534 109 years
 

         b-5)  q0  <  -1

    k = +1  ( spherical geometry )

Formulae:

    D = (c/H0/(-q0)1/2). [ Arccosh (1/b) - Arccosh (1/b/(z+1)) ]      where   b = ((1+q0)/q0)1/2
   D0 = (c/H0/(-1-q0)1/2). [ Arccos b - Arccos b(z+1) ]
   DL = R0 ( 1 + z ) Sin D0/R0          R0 = c/H0/(-1-q0)1/2       R(t) = R0.b Cosh ( H0.t.(-q0)1/2 )      Rmin =  R0.b
   VR = H0 D0
    t0 =  1/H0/(-q0)1/2 Arccosh 1/b  = the time since the minimum of the scale factor
 
 

01  LBL "VAC3"
02  DEG
03  1
04  +
05  STO M
06  X<>Y
07  CHS
08  ENTER^
09  STO N
10  LASTX
11  -
12  /
13  STO O
14  SQRT
15  LASTX
16  1
17  -
18  SQRT
19  +
20  ENTER^
21  R^
22  ST* Y
23  X^2
24  CHS
25  RCL O
26  +
27  SQRT
28  RCL O
29  SQRT
30  +
31  /
32  LN
33  X<> M
34  STO Z
35  RCL O
36  SQRT
37  1/X
38  ST* Y
39  ACOS
40  X<>Y
41  ACOS
42  -
43  STO O
44  SIN
45  R^
46  *
47  X<>Y
48  LN
49  RCL N
50  SQRT
51  ST/ M
52  /
53  RCL N
54  1
55  -
56  SQRT
57  ST/ Z
58  R-D
59  ST/ O
60  RDN
61  SIGN
62  X<> O
63  STO Z
64  137 E8
65  ST* L
66  ST* M
67  ST* Y
68  ST* Z
69  X<> M
70  CLA
71  END

( 107 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:      q0 =  -1.1 ;  z = 2

  -1.1  ENTER^
     2   XEQ "VAC3"   >>>>  D  =  18.458 109 ly
                                  RDN   D0 =  35.699 109 ly
                                 RDN   DL =  95.381 109 ly
                                  RDN   VR =  2.6057
                                LASTX  t0  =  24.408 109 years

***  In these models, there is no Big Bang:    t0 = the time since the minimum of the scale factor R(t)
                                                                    The age of the Universe is actually infinite.

Notes:

-In these Universes, z cannot exceed a maximum redshift:   zmax = (q0/(q0+1))1/2 - 1
-With  q0 = -1.1  ,  zmax = 2.3166

-In fact, there are 2 distances for a given redshift ( one during the contracting phase and another one during the expanding phase )
-This program only gives the distance corresponding to the expanding phase.
 
 

        c) Eddington-Lemaitre Universes ( Critical parameters )
 

-Here,  ( 3L0+2q0-1)3 = 27 L0(q0+L0)2   and  L0 is found by:    6.L0 = ( q0 + 1 )2 - 6.( q0 + 1 ) + 6 + [ ( q0 + 1 )4 -(4/3).( q0 + 1 )3 ]1/2

  Formulae:         with  a = (1/2).(1+q0/L0) -1/3

  H0.D/c = (1/(2a)) sign(-q0).[2a(q0+L0)] -1/2 { 3-1/2 Ln [ ( (z+1+a)1/2+(3a)1/2 ).( (1+a)1/2-(3a)1/2 )/( (z+1+a)1/2-(3a)1/2 )/( (1+a)1/2+(3a)1/2 ) ]
                                                                         + Ln [ ( (z+1+a)1/2-a1/2 ).( (1+a)1/2+a1/2 )2/( (z+1+a)1/2+a1/2 ) ] }

 H0.D0/c =  sign(-q0).[2a(q0+L0)] -1/2  3-1/2 Ln [ ( (z+1+a)1/2+(3a)1/2 ).( (1+a)1/2-(3a)1/2 )/( (z+1+a)1/2-(3a)1/2 )/( (1+a)1/2+(3a)1/2 ) ]

  H0.t0 = (1/(2a)) [2a(q0+L0)] -1/2 { 3-1/2 Ln [ ( (1+a)1/2+(3a)1/2 )/( (1+a)1/2-(3a)1/2 ) ] - Ln ( (1+a)1/2+a1/2 )2 }  if  q0 > 1/2
       t0 = infinity  if  q0 < -1

-All these Universes are spherical  ( k = +1 )
 

Data Registers: /
Flag:  F10
Subroutines: /
 

  01  LBL "LEM"
  02  DEG
  03  CF 10
  04  1
  05  +
  06  STO M
  07  X<>Y
  08  X>0?
  09  SF 10
  10  STO N
  11  LASTX
  12  +
  13  STO Y
  14  4
  15  3
  16  /
  17  -
  18  *
  19  *
  20  *
  21  SQRT
  22  X<> Z
  23  6
  24  -
  25  *
  26  +
  27  6
  28  ST+ Y
  29  /                            here,  X-register = L0    Add  STO 00  after this line if you want to save this number.
  30  STO P                 ( synthetic )
  31  RCL N
  32  X<>Y
  33  ST+ Y
  34  /
  35  PI
  36  INT
  37  1/X
  38  Y^X
  39  ST+ X
  40  1/X                     here,  X-register = a    Add  STO 01  after this line if you want to save this value.
  41  STO Q
  42  RCL M
  43  +
  44  SQRT
  45  STO Y
  46  RCL Q
  47  PI
  48  INT
  49  *
  50  SQRT
  51  STO O
  52  ST+ Z
  53  -
  54  ST/ Y
  55  CLX
  56  SIGN
  57  RCL Q
  58  +
  59  SQRT
  60  RCL X
  61  RCL O
  62  ST+ Z
  63  -
  64  /
  65  ST/ Y
  66  X>0?
  67  LN
  68  STO O
  69  X<>Y
  70  LN
  71  FS? 10
  72  CHS
  73  PI
  74  INT
  75  SQRT
  76  ST/ O
  77  /
  78  RCL M
  79  RCL Q
  80  +
  81  SQRT
  82  RCL X
  83  RCL Q
  84  SQRT
  85  ST- Z
  86  +
  87  /
  88  RCL Q
  89  ENTER^
  90  SIGN
  91  +
  92  SQRT
  93  RCL Q
  94  SQRT
  95  +
  96  X^2
  97  LN
  98  ST- O
  99  X<> L
100  *
101  LN
102  FS? 10
103  CHS
104  +
105  RCL N
106  RCL P
107  +
108  RCL Q
109  ST+ X
110  ST/ O
111  ST/ Z
112  *
113  SQRT
114  ST/ O
115  ST/ Z
116  /
117  RCL P
118  3
119  *
120  RCL N
121  ST+ X
122  +
123  1
124  -
125  SQRT
126  ST/ M
127  R^
128  *
129  R-D
130  SIN
131  ST* M
132  CLX
133  137 E8
134  ST* M
135  ST* O
136  ST* Y
137  ST* Z
138  CLX
139  90
140  TAN
141  FS?C 10
142  X<> O
143  SIGN
144  X<> M
145  X<> Z
146  X<>Y
147  CLA
148  END

( 219 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z             /            DL
           Y             q0            D0
           X              z            D
           L              /            t0

   with  q0 > 1/2   ( Eddington models )
     or   q0 < -1     ( Lemaitre models )

-Such "critical" models don't exist if q0 is between  -1 and 1/2

Example1:      q0 =  -1.1 ;  z = 2      whence   L0 = 1.107976565  ( line 29 )

  -1.1  ENTER^
     2   XEQ "LEM"    >>>>  D  =  17.355 109 ly
                                  RDN   D0 =  32.783 109 ly
                                  RDN   DL =  87.124 109 ly
                                  RDN   VR =  2.3930
                                LASTX  t0  =  9.9999 1099 years = infinity

-If  q0 < -1 ,   z  must be smaller than  2a - 1        ( ymin = 1/(2a) )     here  a = 2.589454154   ( line 40 )
 

Example2:      q0 =  2.375 ;  z = 2      whence   L0 = 1  ( line 29 )

   2.375  ENTER^
     2   XEQ "LEM"    >>>>  D  =   5.0238 109 ly
                                  RDN   D0 =  7.4017 109 ly
                                 RDN   DL =  15.599 109 ly
                                  RDN   VR =  0.5403
                                LASTX  t0  =  5.7825 109 years

-If  q0 > 1/2  ,  ymax = 1/(2a)  ,   here   a = 1/3

-Writing  SIZE 000 programs is always an interesting challenge,
  but several bytes may be saved if you use standard registers ...
 

        d) Euclidean Universes
 

-We now assume k = 0   i-e   3.L0 + 2.q0 - 1 = 0

   q0 = 1/2  is Einstein-de Sitter model
   q0 = -1   is de Sitter model

-For  other q0-values,  the comoving distance requires elliptic integrals
  but the light-travel time distance D and the age of the Universe may be expressed in terms of elementary functions:
 

           d-1)  -1 < q0 < 1/2
 

Formulae:    H0.t0 =  (3-6q0) -1/2  Ln ( (1+a)1/2 + a1/2 )2         with  a = (1-2q0)/(2+2q0)

    H0.D/c = (3-6q0) -1/2  Ln [ ( ( (z+1)3 + a )1/2 - a1/2 ) ( ( (z+1)3 + a )1/2 + a1/2 ) -1 ( (1+a)1/2 + a1/2 )2 ]
 

01  LBL "EUCL"
02  X<>Y
03  ENTER^
04  ST+ X
05  1
06  ST+ Z
07  ST+ T
08  -
09  CHS
10  STO M
11  X<>Y
12  ST+ X
13  /
14  ENTER^
15  R^
16  3
17  Y^X
18  +
19  SQRT
20  ENTER^
21  R^
22  SQRT
23  ST- Z
24  +
25  /
26  1
27  R^
28  ST+ Y
29  SQRT
30  X<>Y
31  SQRT
32  +
33  X^2
34  ST* Y
35  LN
36  X<>Y
37  LN
38  0
39  X<> M
40  3
41  *
42  SQRT
43  137 E8
44  X<>Y
45  /
46  ST* Z
47  *
48  END

( 71 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Y            q0            t0
           X              z            D

Example:     q0 = -0.7  &  z = 7      ( whence  L0 = 0.8 )

    -0.7  ENTER^
       7   XEQ "EUCL"   >>>>   D = 1.3840  1010 ly
                                    X<>Y   t0 = 1.4742  1010 years
 

           d-2)  q0 > 1/2
 

Formulae:    H0.t0 =  (-3+6q0) -1/2 [ PI - 2.Arctan ( (1+a)1/2  (-a) -1/2 ) ]      with  a = (1-2q0)/(2+2q0)

    H0.D/c = 2.(-3+6q0) -1/2   [ Arctan ( ( (z+1)3 + a )1/2 (-a) -1/2 ) - Arctan ( ( 1 + a )1/2 (-a) -1/2  ]
 
 

01  LBL "EUCL2"
02  DEG
03  X<>Y
04  ENTER^
05  ST+ X
06  1
07  ST+ Z
08  ST+ T
09  -
10  STO M
11  X<>Y
12  ST+ X
13  /
14  ENTER^
15  R^
16  3
17  Y^X
18  -
19  CHS
20  1
21  R^
22  ST- Y
23  ST/ Z
24  /
25  SQRT
26  ATAN
27  X<>Y
28  SQRT
29  ATAN
30  X<>Y
31  ST- Y
32  D-R
33  CHS
34  ST+ X
35  PI
36  +
37  X<>Y
38  ST+ X
39  D-R
40  0
41  X<> M
42  3
43  *
44  SQRT
45  137 E8
46  X<>Y
47  /
48  ST* Z
49  *
50  END

( 76 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             q0            t0
           X              z            D

Example:     q0 = 2  &  z = 7      ( whence  L0 = -1 )

       2   ENTER^
       7   XEQ "EUCL"   >>>>   D = 6.8878  109 ly
                                    X<>Y   t0 = 7.1733  109 years
 
 

        e) More General Models
 

           e-1)  Numerical Integration ( program#1 )
 

-Three integrals are needed to compute  t0 , D and D0
-The following program uses the change of variable  y = u2 and Gauss-Legendre 3-point formula.
 

Data Registers:   R00 thru R16        R00 thru R08 are used by  "GL3"

             R09:  z + 1               R12:   2(q0 + L0)                 R15-R16: temp
             R10:  q0                   R13:   2q0 + 3L0 - 1
             R11:  L0                   R14:   (1/H0).t0

 when the program stops:   R00 = k     ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
                                             R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )

Flags:  F00  set flag F00  if you don't want to re-calculate the age of the Universe t0
            F10
Subroutine:   "GL3"  Gauss-Legendre 3-point formula  ( cf "Numerical Integration for the HP-41" )
 
 

  01  LBL "COSMO"
  02  DEG
  03  1
  04  STO 02
  05  +
  06  STO 09
  07  "S"
  08  ASTO 00
  09  SF 10
  10  FS? 00
  11  GTO 01
  12  CLX
  13  STO 01
  14  +
  15  STO 10
  16  X<>Y
  17  STO 11
  18  +
  19  X<0?
  20  LN                    you'll get   DATA ERROR  if the density of matter is negative  i-e  q0 + L0 < 0
  21  STO 12
  22  ST+ X
  23  +
  24  1
  25  -
  26  STO 13
  27  3
  28  Y^X
  29  RCL 12
  30  ST+ 12
  31  X^2
  32  27
  33  *
  34  RCL 11
  35  *
  36  X>Y?
  37  GTO 00
  38  RCL 10
  39  1
  40  +
  41  X>0?
  42  GTO 00
  43  90
  44  TAN                         if   ( 3L+2q-1)3 >= 27 L(q+L)2   and  q <= -1  there is no Big Bang  and  t0 = infinity ( 9.999999999 E99 for the HP-41 )
  45  STO 14
  46  GTO 01
  47  LBL 00
  48  XEQ 02
  49  STO 14
  50  LBL 01
  51  RCL 09
  52  SQRT
  53  1/X
  54  STO 01
  55  XEQ 02
  56  STO 15
  57  CF 10
  58  XEQ 02
  59  STO 02
  60  STO 04
  61  RCL 13
  62  ABS
  63  SQRT
  64  X=0?
  65  SIGN
  66  1/X
  67  STO 06
  68  /
  69  ENTER^
  70  E^X-1
  71  LASTX
  72  CHS
  73  E^X-1
  74  -
  75  2
  76  /
  77  RCL Y
  78  R-D
  79  SIN
  80  RCL 13
  81  X#0?
  82  ENTER^
  83  X>0?
  84  ENTER^
  85  GTO 04
  86  LBL 02                 Lines 86 to 103  use Gauss Legendre 3-point formula
  87   E99                      the number n of subintervals ( in register R03 ) is doubled until 2 consecutive rounded values are equal
  88  STO 16                 thus, the accuracy is controlled by the display format
  89  SIGN
  90  STO 03                Actually, for our Universe ( q0 of the order of  -0.5 and L0 of the order of 0.6 ),
  91  LBL 03                 the results have at least 6 significant digits if n = 4 and lines 87 to 103 may be replaced by  4  STO 03  XEQ "GL3"  RTN
  92  RCL 03
  93  ST+ 03                 With our Universe, n = 10 produces almost 10 significant figures!
  94  XEQ "GL3"
  95  VIEW X
  96  RND
  97  LASTX
  98  X<> 16
  99  RND
100  X#Y?
101  GTO 03
102  RCL 04
103  RTN
104  LBL "S"
105  X^2
106  ENTER^
107  ENTER^
108  X^2
109  RCL 11
110  *
111  RCL 13
112  -
113  *
114  RCL 12
115  +
116  SQRT
117  1/X
118  ST+ X
119  FS? 10
120  *
121  RTN
122  LBL 04
123  R^
124  RCL 06
125  *
126  RCL 09
127  *
128  STO 03
129  RCL 13
130  X#0?
131  SIGN
132  STO 00
133  RCL 14
134  STO 05
135  RCL 15
136  STO 01
137  137 E8
138  ST* 01
139  ST* 02
140  ST* 03
141  SF 25
142  ST* 05
143  ST* 06
144  RCL 05
145  SIGN
146  RCL 04
147  RCL 03
148  RCL 02
149  RCL 01
150  CLD
151  CLA
152  CF 25
153  END

( 212 bytes / SIZE 017 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z             L0            DL
           Y            q0            D0
           X              z            D
           L              /            t0

Example1:     L0 = 0.7  ;  q0 =  -0.6   ;   z = 7

*** These values correspond to  (Omega)vac = L0 = 0.7  &  (Omega)mat =  2(L0+q0) = 0.2   whence (rho)mat = 3H02(L0+q0)/(4PI.G) ~ 1.91 10-27 kg/m3

     0.7  ENTER^
   -0.6  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.3283 1010 ly         = R01            ( in 93s in FIX 4 )             The successive approximations
                                      RDN   D0 =  3.0246 1010 ly         = R02                                                     of the 3 integrals (  for  t0  D0  DL )
                                      RDN   DL =  2.6210 1011 ly         = R03                                                     are displayed during the calculations.
                                      RDN   VR =  2.2077                     = R04
                                    LASTX  t0  =  1.4168 1010 years    = R05

-We also have    R06 = 4.3323 1010 = current scale factor = R(t0) = R0
            and        R00 = k = -1 , meaning the Universe is hyperbolic.

-Registers R10 thru R14 will be initialized after executing "COSMO" once.
-Set flag F00 and simply place z in X-register to compute other distances
 with the same cosmological parameters and the same accuracy.

-This program may be used to draw a graphic representation of R(t):
 (  add FS? 04  RTN  after line 56  and set flags  F00  and  F04  in order to calculate H0.( t0 - t )  only )

 for example:    z = 7  means  y = 1/8  (  z = 1/y - 1 ) and  light-travel time = 1.3283 1010 =  t0 - t
  whence  t = 8.85 1010  or  H0.t = 0.0646  ( if  y = 1/8 )

-More directly, add   LBL "IG"  after line 86  LBL 02  and load this short routine:

01  LBL "Y-H0T"                      In this program, the origine of time ( t = 0 ) is the Big Bang.
02  X=0?                                   Of course, it will not work if there is no Big Bang!
03  RTN                                    In this case, replace lines 02 to 08 by
04  SF 10
05  SQRT                                  SF 10   SQRT   STO 02   1   STO 01
06  STO 02
07  CLX                                    Thus, the origin of time will be "now"
08  STO 01
09  "S"
10  ASTO 00
11  XEQ "IG"
12  CLA
13  END

-Initialize registers R11 to R13 ( this is done if you have already executed "COSMO" )
 and with our example:    q0 =  -0.6  ;  L0 = 0.7    it yields  ( in FIX 3 )

  For y = 3    3  XEQ "Y-H0T"  >>>>  2.267  and similarly:
 
      y       0     0.1    0.25      yI    0.75      1    1.25     1.5    1.75      2    2.5      3     10
   H0.t       0   0.046   0.178   0.495   0.765   1.034   1.266   1.466   1.640   1.793   2.053   2.267   3.700

 

  y = R(t)/R0
        |
    3  |                                                                   *
        |
        |                                                               *
        |
    2  |                                                        *
        |                                                    *
        |                                               *
        |                                         *
    1  |                                  * now  ( y = 1 )
        |                         *
        |                * Infl
        |     *
        |*----------|----------|----------|----------|-------  H0.t
       0              0.5             1            1.5            2

 Here,  yI = ( 1 + q0/L0)1/3 = 0.522757959      I  is a point of inflexion:     q > 0  for  y < yI   and  q < 0  for y > yI
 

Example2:     L0 = 2/3  ;  q0 =  -0.5   ;   z  = 7

*** These values correspond to  (Omega)vac = L0 = 2/3  &  (Omega)mat =  2(L0+q0) = 1/3   whence (rho)mat = 3H02(L0+q0)/(4PI.G) ~ 3.19 10-27 kg/m3

      2    ENTER^  3  /
   -0.5  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.2123 1010 ly         = R01            ( in 97s in FIX 4 )
                                      RDN   D0 =  2.6605 1010 ly         = R02
                                      RDN   DL =  2.1284 1011 ly         = R03
                                      RDN   VR =  1.9420                     = R04
                                    LASTX  t0  =  1.2822 1010 years    = R05

-We also have    R06 = 1.3700 1010 = current scale factor = R(t0) = R0
            and        R00 = k = 0 , we have an Euclidean Universe.
 

Example3:     L0 = 0.6  ;  q0 =  -0.3   ;   z  = 7

*** These values correspond to  (Omega)vac = L0 = 0.4  &  (Omega)mat =  2(L0+q0) = 0.6   whence (rho)mat = 3H02(L0+q0)/(4PI.G) ~ 5.74 10-27 kg/m3

     0.6  ENTER^
   -0.3  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.0689 1010 ly         = R01            ( in 120s in FIX 4 )
                                      RDN   D0 =  2.2467 1010 ly         = R02
                                      RDN   DL =  1.6405 1011 ly         = R03
                                      RDN   VR =  1.6399                     = R04
                                    LASTX  t0  =  1.1217 1010 years    = R05

-We also have    R06 = 3.0634 1010 = current scale factor = R(t0) = R0
            and        R00 = k = +1 , the Universe is spherical.
 

Notes:

-If there is no Big Bang,  the minimum  ymin = Rmin/R0 may be found by solving the cubic equation:   L0.y3 + ( 1-2q0-3L0 ).y +  ( 2q0+2.L0 )  = 0

  For example, if    q0 = -1.1 & L0 = 1.102   the cubic  1.102 y3 - 0.106 y + 0.004 = 0  has 3 real roots:  0.289201975 ; 0.038320885 ; -0.327522859
  and   ymin = 0.289201975  which corresponds to  zmax = 1/ymin - 1 = 2.457791053

-If  z is close to zmax , "COSMO"  will be very slow:  thus, in FIX 3

  1.102  ENTER^
   -1.1   ENTER^
    2.45  XEQ "COSMO"  yields  D =  2.459 1010                 in  1h16mn      for only 2 integrals !
                                                   D0 = 5.599 1010
                                                   DL = 1.410 1011
                                                   VR = 4.087           R0 = 4.208 1010

-For similar reasons, long execution time is to be expected if you are looking for the half-period of an oscillating universe ( Domains 1&2 )
-"Plateau" models ( when the parameters L and q are close to their "critical" values ) are a third difficult case.

-In all these 3 cases, the integrand is ( or is almost ) singular and a good 41-emulator in turbo mode will be useful!
 

            e-2)  Numerical Integration ( program#2 )
 

-The program above gives accurate results in most cases, especially for our Universe.
-If, however, you want to explore other models, other changes of variable may be better. "COSMO2" uses the following ones:

   a) If  y  has a positive minimum ( bounce models )  we set:   y = ymin cosh2 v
   b) If  y  has a finite maximum ( oscillating models )  we set:   y = ymax sin2 v
   c) For "plateau" models, we set   y = b + e.sinh v  where  b+i.e , b-i.e  are the complex roots of the radicand ( e > 0 ).

-This last change of argument is put forward with little enthusiasm because numerical integration remains difficult for plateau models.
-However, it performs much better than "COSMO", and even better if you are using the integrator listed in the PPC ROM user manual.
 

Data Registers:   R00 thru R22

      R10:  z + 1             R13:  2(q0 + L0)                   R16:   ymax                               R19:  proportional to D
      R11:  q0                 R14:   -2q0 - 3L0 + 1            R17:   proportional to t0            R20-R21-R22: temp
      R12:  L0                 R15:  ymin                             R18:   proportional to T

 when the program stops:   R00 = k       ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
                                             R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )
                                             R07 = Rmin  ,  R08 = Rmax  ,   R09 = T = Period of an oscillating Universe ( 0 otherwise )

-Here, for bounce Universes,  t0 is reckoned from the minimum of the scale factor R(t)  ( the age of the Universe is infinite )

Flags:  F00  set flag F00  if you don't want to (re-)calculate the age of the Universe t0and the period T
            F01  is used by "P3" to indicate complex roots
            F10

Subroutines:   "P3"    Cubic equations   ( cf  "Polynomials for the HP-41" )
                         "GL3"  Gauss-Legendre 3-point formula  ( cf "Numerical Integration for the HP-41" )
 

  01  LBL "COSMO2"
  02  DEG
  03  CF 01
  04  SF 10
  05  1
  06  STO 02
  07  +
  08  STO 10
  09  CLX
  10  STO 01
  11  STO 18
  12  +
  13  STO 11
  14  X<>Y
  15  STO 12
  16  +
  17  ST+ X
  18  X<0?
  19  LN
  20  STO 13
  21  1
  22  RCL Y
  23  -
  24  R^
  25  -
  26  STO 14
  27  R^
  28  X=0?
  29  GTO 00
  30  CLX
  31  X<> Z
  32  XEQ "P3"
  33  GTO 01
  34  LBL 00
  35  RDN
  36  CHS
  37  X<=0?
  38  CLST
  39  X#0?
  40  /
  41  LBL 01
  42  STO 03
  43  RDN
  44  STO 21
  45  FS? 01
  46  CLX
  47  STO 04
  48  X<>Y
  49  STO 22
  50  FS? 01
  51  CLX
  52  STO 05
  53  1
  54  RCL 03
  55  X<Y?
  56  GTO 00
  57  CLX
  58  RCL 04
  59  X<Y?
  60  GTO 00
  61  CLX
  62  RCL 05
  63  X>Y?
  64  CLX
  65  LBL 00
  66  X<0?
  67  CLX
  68  STO 15
  69  CLX
  70  RCL 05
  71  X>Y?
  72  GTO 00
  73  CLX
  74  RCL 04
  75  X>Y?
  76  GTO 00
  77  CLX
  78  RCL 03
  79  X>Y?
  80  GTO 00
  81  90
  82  TAN
  83  LBL 00
  84  STO 16
  85  90
  86  TAN
  87  -
  88  RCL 15
  89  +
  90  X=0?
  91  GTO 03
  92  LASTX
  93  X=0?
  94  GTO 02
  95  "S1"
  96  ASTO 00
  97  SQRT
  98  1/X
  99  RCL 15
100  1/X
101  1
102  -
103  SQRT
104  +
105  LN
106  STO 02
107  RCL 03
108  RCL 04
109  X#Y?
110  GTO 00
111  90
112  TAN
113  GTO 01
114  LBL 00
115  FC? 00
116  XEQ 05
117  LBL 01
118  FC? 00
119  STO 17
120  RCL 10
121  RCL 15
122  *
123  1/X
124  SQRT
125  LASTX
126  1
127  -
128  SQRT
129  +
130  LN
131  GTO 04                      ( theoretically a three-byte GTO )
132  LBL 02
133  "S2"
134  ASTO 00
135  90
136  STO 02
137  RCL 03
138  RCL 04
139  X=Y?
140  GTO 01
141  FC? 00
142  XEQ 05
143  ST+ X
144  FC? 00
145  STO 18
146  LBL 01
147  RCL 16
148  SQRT
149  1/X
150  ASIN
151  STO 02
152  FC? 00
153  XEQ 05
154  FC? 00
155  STO 17
156  RCL 10
157  RCL 16
158  *
159  SQRT
160  1/X
161  ASIN
162  GTO 04
163  LBL 03
164  RCL 11
165  1
166  +
167  FS? 01
168  X>0?
169  GTO 03
170  20
171  1/X
172  RCL 22
173  ABS
174  STO 22
175  X>Y?
176  GTO 03
177  "S3"
178  ASTO 00
179  RCL 21
180  X<>Y
181  /
182  ENTER^
183  X^2
184  1
185  +
186  SQRT
187  +
188  LN
189  CHS
190  STO 01
191  1
192  RCL 21
193  -
194  RCL 22
195  /
196  ENTER^
197  X^2
198  1
199  +
200  SQRT
201  +
202  LN
203  STO 02
204  FC? 00
205  XEQ 05
206  FC? 00
207  STO 17
208  RCL 10
209  1/X
210  RCL 21
211  -
212  RCL 22
213  /
214  ENTER^
215  X^2
216  1
217  +
218  SQRT
219  +
220  LN
221  GTO 04
222  LBL 03
223  "S"
224  ASTO 00
225  RCL 11
226  1
227  +
228  RCL 12
229  X=Y?
230  X#0?
231  GTO 00
232  90
233  TAN
234  GTO 01
235  LBL 00
236  FC? 00
237  XEQ 05
238  LBL 01
239  FC? 00
240  STO 17
241  RCL 10
242  SQRT
243  1/X
244  LBL 04
245  STO 01
246  XEQ 05
247  STO 19
248  CF 10
249  XEQ 05
250  STO 02
251  STO 20
252  RCL 14
253  ABS
254  SQRT
255  X=0?
256  SIGN
257  1/X
258  STO 03
259  STO 06
260  STO 07
261  STO 08
262  /
263  ENTER^
264  E^X-1
265  LASTX
266  CHS
267  E^X-1
268  -
269  2
270  /
271  RCL Y
272  R-D
273  SIN
274  RCL 14
275  X#0?
276  ENTER^
277  X<0?
278  ENTER^
279  GTO 07                      ( theoretically a three-byte GTO )
280  LBL 05
281   E99
282  STO 20
283  SIGN
284  STO 03
285  LBL 06
286  RCL 03
287  ST+ 03
288  XEQ "GL3"
289  VIEW X
290  RND
291  LASTX
292  X<> 20
293  RND
294  X#Y?
295  GTO 06
296  RCL 04
297  RTN
298  LBL "S"
299  X^2
300  ENTER^
301  ENTER^
302  X^2
303  RCL 12
304  *
305  RCL 14
306  +
307  *
308  RCL 13
309  +
310  SQRT                        5 bytes may be saved if you add    LBL 00   after line  336
311  1/X                            and if you replace lines 310 to 315 by  GTO 00
312  ST+ X
313  FS? 10
314  *
315  RTN
316  LBL "S1"
317  E^X
318  ENTER^
319  1/X
320  +
321  2
322  /
323  X^2
324  RCL 15
325  *
326  ENTER^
327  STO Z
328  RCL 15
329  +
330  *
331  RCL 12
332  *
333  RCL 13
334  RCL 15
335  /
336  -
337  SQRT
338  1/X
339  ST+ X
340  FS? 10
341  *
342  RTN
343  LBL "S2"
344  SIN
345  X^2
346  RCL 16
347  *
348  ENTER^
349  STO Z
350  RCL 16
351  +
352  *
353  RCL 12
354  CHS
355  *
356  RCL 13
357  RCL 16
358  /
359  +
360  SQRT
361  1/X
362  ST+ X
363  D-R
364  FS? 10
365  *
366  RTN
367  LBL "S3"
368  E^X
369  ENTER^
370  1/X
371  -
372  2
373  /
374  RCL 22
375  *
376  RCL 21
377  +
378  STO Y
379  LASTX
380  ST+ X
381  +
382  RCL 12
383  *
384  FS? 10
385  GTO 00
386  *
387  SQRT
388  1/X
389  RTN
390  LBL 00
391  /
392  SQRT
393  RTN
394  LBL 07
395  R^
396  RCL 10
397  *
398  ST* 03
399  RCL 14
400  X#0?
401  SIGN
402  CHS
403  STO 00
404  RCL 17
405  STO 05
406  RCL 19
407  STO 01
408  RCL 15
409  ST* 07
410  SF 24
411  RCL 16
412  ST* 08
413  RCL 18
414  STO 09
415  9
416  137 E8
417  LBL 08
418  ST* IND Y
419  DSE Y
420  GTO 08
421  RCL 05
422  SIGN
423  RCL 20
424  STO 04
425  RCL 03
426  RCL 02
427  RCL 01
428  CLA
429  CLD
430  CF 24
431  END

( 602 bytes / SIZE 023 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z             L0            DL
           Y             q0            D0
           X              z            D
           L             /            t0

 

Example1:     L0 = 1.102  ;  q0 =  -1.1   ;   z  = 2.45

    1.102    ENTER^
    -1.1      ENTER^
     2.45    XEQ "COSMO2"  >>>>  D  =  2.4589 1010 ly           = R01            ( in 110s in FIX 3 instead of more than 1 hour with "COSMO" )
                                            RDN   D0 =  5.5990 1010 ly           = R02
                                             RDN   DL =  1.4101 1011 ly           = R03
                                             RDN   VR =  4.0869                       = R04
                                             LASTX  t0  =  2.5501 1010 years    = R05

-We also have    R00 = k = 1   ( Spherical Universe )
                          R06 = 4.2079 1010 = current scale factor = R(t0) = R0
                          R07 = 1.2169 1010 = Rmin  ( bounce model )
                          R08 = 9.9999 1099 = Rmax = infinity
                          R09 = 0  ( non-oscillating Universe )
 

Example2:     L0 = 0.2  ;  q0 =  2   ;   z  = 7

     0.2  ENTER^
       2   ENTER^
       7     R/S  >>>>  D  =  6.1763 109 ly          = R01            ( in 158s in FIX 3 )
                     RDN   D0 =  1.1392 1010 ly         = R02
                     RDN   DL =  5.7763 1010 ly         = R03
                     RDN   VR =  0.8315                     = R04
                   LASTX  t0  =  6.3750 109 years     = R05

-We also have    R00 = k = +1  ( Spherical Universe )
                          R06 = 7.2205 109 = current scale factor = R(t0) = R0
                          R07 = 0 = Rmin
                          R08 = 9.8405 109 = Rmax
                          R09 = 3.5658 1010 = T = time between a Big Bang and a Big Crunch.
 

Example3:     L0 = 1.108  ;  q0 =  -1.1   ;   z  = 7

    1.108    ENTER^
    -1.1      ENTER^
       7          R/S   >>>>  D  =  6.7540 1010 ly           = R01            ( in 11mn25s in FIX 3  )
                          RDN   D0 =  2.8508 1011 ly           = R02
                          RDN   DL =  2.6910 1011 ly           = R03
                          RDN   VR =  20.8088                     = R04
                          LASTX  t0  =  7.2650 1010 years    = R05

        and            R00 = k = 1   ( Spherical Universe )
                          R06 = 3.8905 1010 = current scale factor = R(t0) = R0
                          R07 = 0  = Rmin
                          R08 = 9.9999 1099 = Rmax = infinity
                          R09 = 0  ( non-oscillating Universe )

-This is a "plateau model" :  it spends a long time in a quasistatic state.
-  D and D0 are computed ( relatively ) fast ( provided z is not too large ), but t0 requires much more iterations.
-For  q0 =  -1.1  the critical  L0 is approximately  1.107976565  ( cf Lemaitre models )
 
 

            e-3)  Elliptic Integrals ( Big Bang Universes only )
 
 

-  t0 , D and D0 may be computed by Legendre elliptic integrals ( cf the last 2 references below )
-  "COSMO3"  uses Carlson elliptic integrals instead: the formulas are simpler.

              D  = (c/H0)    §1z+1 du / u / ( 2(q0+L0 ).u3 + (1-2q0-3L0 ).u2 + L0 )1/2         u = 1/y      ;     t0 corresponds to z = infinity
              D0 = (c/H0)    §1z+1 du / ( 2(q0+L0 ).u3 + (1-2q0-3L0 ).u2 + L0 )1/2

-If the cubic    2(q0+L0 ).u3 + (1-2q0-3L0 ).u2 + L0 = 0   has 3 real roots  a , b , c  we have:

       D = (c/3/H0) (2/(q0 + L0 ))1/2  [  R J ( -a+1 , -b+1 , -c+1 , 1 )  -  R J ( -a+z+1 , -b+z+1 , -c+z+1 , z+1 )  ]
       D0 = (c/H0) (2/(q0 + L0 ))1/2   [  R F ( -a+1 , -b+1 , -c+1 )  -  R F ( -a+z+1 , -b+z+1 , -c+z+1 )  ]
 

-If the cubic    2(q0+L0 ).u3 + (1-2q0-3L0 ).u2 + L0 = 0   has only 1 real root a  and a pair of complex conjugate roots  b + i.c , b - i.c  we have:

      D = (c/3/H0) (2/(q0 + L0 ))1/2  [  R J ( -a+1 , -b+1+i.c , -b+1-i.c , 1 )  -  R J ( -a+z+1 , -b+z+1+i.c , -b+z+1-i.c , z+1 )  ]
      D0 = (c/H0) (2/(q0 + L0 ))1/2   [  R F ( -a+1 , -b+1+i.c , -b+1-i.c )  -  R F ( -a+z+1 , -b+z+1+i.c , -b+z+1-i.c )  ]

-This program doesn't work if there is no big bang or if the Universe is empty.
 

Data Registers:   R00 thru R21

            R13:  -a                   R16:   z + 1                         R19:  proportional to t0
             R14:  -b                   R17:   q0 + L0                     R20:  ymax
             R15:  -c                   R18:   2q0 + 3L0 - 1            R21:  proportional to T

 when the program stops:   R00 = k         ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
                                             R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )
                                             R07 = Rmax  ,  R08 = T = Period of an oscillating Universe ( 0 otherwise )

Flags:  F00  set flag F00  if you don't want to re-calculate the age of the Universe t0and the period T
            F01  is used by "P3" to indicate complex roots

Subroutines:   "P3"              ( cf  "Polynomials for the HP-41" )
                         "RF" & "RJ"  Carlson elliptic integrals of the first and third kind       ( cf "Jacobian Elliptic Functions and Elliptic Integrals for the "HP-41" )
                         "RFZ" & "RJZ"  ---------------------------------------------- complex arguments ----------------------------------------------------
 

  01  LBL "COSMO3"
  02  1
  03  +
  04  STO 16
  05  FS? 00
  06  GTO 04                      ( theoretically a three-byte GTO )
  07  CLX
  08  +
  09  +                                 add   X<0?   LN    after this line if you want to get "DATA ERROR" for a negative density of matter
  10  STO 17
  11  ST+ X
  12  STO Z
  13  +
  14  1
  15  -
  16  STO 18
  17  0
  18  STO 21
  19  R^
  20  CHS
  21  XEQ "P3"
  22  STO 13
  23  RDN
  24  STO 14
  25  X<>Y
  26  STO 15
  27  X<>Y
  28  R^
  29  1
  30  XEQ 02
  31  STO 19
  32  90                                lines 32 to 69 may be deleted if you don't want to calculate Rmax and  T
  33  TAN
  34  STO 20
  35  SIGN
  36  CHS
  37  FS? 01
  38  GTO 00
  39  RCL 15
  40  X>Y?
  41  GTO 01
  42  CLX
  43  RCL 14
  44  X>Y?
  45  GTO 01
  46  X<>Y
  47  LBL 00
  48  RCL 13
  49  X<Y?
  50  CLX
  51  LBL 01
  52  CHS
  53  X<=0?
  54  GTO 04
  55  1/X
  56  STO 20
  57  RCL 15
  58  RCL 14
  59  X<0?                       lines 59 to 64 are only useful if the Universe is an Eddington model ( finite Rmax but no Big Crunch )
  60  X#Y?
  61  FS? 30
  62  FS? 01
  63  FS? 30
  64  GTO 04
  65  RCL 13
  66  LASTX
  67  XEQ 02
  68  ST+ X
  69  STO 21
  70  GTO 04
  71  LBL 02
  72  FC? 01
  73  ST+ T
  74  ST+ Z
  75  ST+ Y
  76  RDN
  77  FS? 01
  78  XEQ "RJZ"
  79  FC? 01
  80  XEQ "RJ"
  81  RTN
  82  LBL 03
  83  FC? 01
  84  ST+ T
  85  ST+ Z
  86  +
  87  FS? 01
  88  XEQ "RFZ"
  89  FC? 01
  90  XEQ "RF"
  91  RTN
  92  LBL 04
  93  RCL 15
  94  RCL 14
  95  RCL 13
  96  RCL 16
  97  XEQ 02
  98  RCL 19
  99  X<>Y
100  -
101  STO 10
102  RCL 15
103  RCL 14
104  RCL 13
105  RCL 16
106  XEQ 03
107  STO 12
108  RCL 15
109  RCL 14
110  RCL 13
111  1
112  XEQ 03
113  RCL 12
114  -
115  STO 02
116  RCL 10
117  STO 01
118  RCL 19
119  STO 05
120  RCL 20
121  STO 07
122  RCL 21
123  STO 08
124  2
125  RCL 17
126  /
127  SQRT
128  ST* 02
129  3
130  /
131  ST* 01
132  ST* 05
133  ST* 08
134  RCL 18
135  X#0?
136  SIGN
137  STO 00
138  RCL 02
139  STO 09
140  RCL 18
141  ABS
142  SQRT
143  X=0?
144  SIGN
145  1/X
146  STO 06
147  SF 24
148  ST* 07
149  /
150  ENTER^
151  E^X-1
152  LASTX
153  CHS
154  E^X-1
155  -
156  2
157  /
158  RCL Y
159  R-D
160  SIN
161  RCL 00
162  X#0?
163  ENTER^
164  X>0?
165  ENTER^
166  R^
167  RCL 06
168  *
169  RCL 16
170  *
171  STO 03
172  8
173  137 E8
174  LBL 05
175  ST* IND Y
176  DSE Y
177  GTO 05
178  RCL 05
179  SIGN
180  RCL 09
181  STO 04
182  RCL 03
183  RCL 02
184  RCL 01
185  CF 24
186  END

( 286 bytes / SIZE 022 )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z             L0            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example1:     L0 = 0.7  ;  q0 =  -0.6   ;   z  = 7

     0.7  ENTER^
   -0.6  ENTER^
      7   XEQ "COSMO3"  >>>>  D  =  1.328255073 1010 ly         = R01            ( in 140s )      flag F01 is set
                                        RDN   D0 =  3.024557313 1010 ly         = R02
                                        RDN   DL =  2.621046346 1011 ly         = R03
                                        RDN   VR =  2.207706068                     = R04
                                      LASTX  t0  =  1.416778742 1010 years    = R05

-We also have    R00 = k = -1 , meaning the Universe is hyperbolic.
                          R06 = 4.332320394 1010 = current scale factor = R(t0) = R0
                          R07 = 9.999999999 1099 = infinity = Rmax
                          R08 = 0  ( non oscillating Universe )
 

Example2:     L0 = 0.2  ;  q0 =  2   ;   z  = 7

     0.2  ENTER^
       2   ENTER^
       7   XEQ "COSMO3"  >>>>  D  =  6.176349709 109 ly          = R01            ( in 145s )      flag F01 is clear
                                        RDN   D0 =  1.139222269 1010 ly         = R02
                                        RDN   DL =  5.776287398 1010 ly         = R03
                                        RDN   VR =  0.831549101                     = R04
                                      LASTX  t0  =  6.375024998 109 years     = R05

-We also have    R00 = k = +1  ( Spherical Universe )
                          R06 = 7.220533991 109 = current scale factor = R(t0) = R0
                          R07 = 9.840506634 109 = Rmax
                          R08 = 3.565763808 1010 = T = time between a Big Bang and a Big Crunch.

-Registers R13 thru R21 will be initialized after executing "COSMO3" once.
-Set flag F00 , do not change flag F01 and simply place z in X-register to compute other distances with the same cosmological parameters.
 

Note:      With  L0 = 1.108  &  q0 =  -1.1  "COSMO3"  yields  t0  = 7.264885687 1010 years
 

        f) Non-Vanishing Pressure
 

-We now assume that the Universe is made of dust & radiation.
-In this case, Einstein's equations yields:

          2.q + 2.L = 8.PI.G.(rho)/(3H2) + ¶                 where      ¶ = 8.PI.G.prad/(c2H2)  is the pressure parameter = (Omega)rad
     2.q + 3.L - 1 = k.c2/(H2R2) + ¶                             and        prad is the pressure of radiation

   (rho) = (rho)mat + 3.prad/c2   and according to Stefan's law:   3prad = (rho)rad c2 = a.T4
                                                                                   where          a = Stefan's Constant  and  T = temperature of the Cosmic Microwave Background

-If we assume there is no creation/annihilation of matter,  (rho)mat R3 = constant  and   prad R4 = constant.
-The integrals become:

     D  = (c/H0)  §y(em)1  y. [  L0.y4 + ( 1+¶0-2q0-3L0 ).y2 + 2.( q0+L0- ¶0 ).y + ¶0  ] -1/2 dy           y(em) = y at the instant of emission     z + 1 = R0/R = 1/yem
     D0 = (c/H0)  §y(em)1  [  L0.y4 + ( 1+¶0-2q0-3L0 ).y2 + 2.( q0+ L0- ¶0 ).y + ¶0  ] -1/2 dy              D0 = comoving radial distance

-We now have  (Omega)mat = 2(q0+L0-¶0)
 

            f-1)  Pure Radiation
 

-Then, (Omega)mat = 0 whence  ¶0 = q0+L0
-Such models have been studied by Tolman.
-There are many cases, and the following program only works for Big Bang Universes with a positive Cosmological Constant.
-The comoving radial distance involves elliptic integrals and "TOL" calculates the light-travel time Distance D,
  the age of the Universe t0 , the current radius of the Universe R0 and the constant  k
  ( D & t0 may be expressed in terms of elementary functions )

Formula:

  H0.(t0 - te ) = 1/(2.L01/2) Ln { [ 1 + (1-q0-2.L0)/(2.L0) + 1/L01/2 ] / [ 1/(z+1)2 + (1-q0-2.L0)/(2.L0) + ( 1/(z+1)4 + ( (1-q0-2.L0)/(z+1)2 + q0+ L0 )/L0 )1/2 ] }

  ( replace z by +infinity to get  t0 )

Data Registers: /
Flags: /
Subroutines: /
 

01  LBL "TOL"
02  1
03  +
04  X^2
05  1/X
06  STO N
07  CLX
08  +
09  +
10  STO M
11  +
12  STO O
13  SIGN
14  ST- O
15  X<>Y
16  SQRT
17  1/X
18  +
19  RCL O
20  RCL N
21  *
22  RCL M
23  X<>Y
24  -
25  R^
26  /
27  RCL N
28  X^2
29  +
30  SQRT
31  RCL O
32  R^
33  ST/ M
34  STO P           ( synthetic )
35  ST+ X
36  /
37  ST- Y
38  ST- Z
39  ST- T
40  X<> N
41  +
42  /
43  LN
44  X<> M
45  SQRT
46  RCL N
47  -
48  /
49  LN
50  RCL P
51  SQRT
52  ST+ X
53  ST/ M
54  /
55  RCL O
56  ABS
57  LASTX
58  X#0?
59  SIGN
60  RDN
61  X=0?
62  SIGN
63  SQRT
64  1/X
65  STO Z
66  CLX
67  137 E8
68  ST* M
69  ST* Z
70  *
71  RCL M
72  CLA
73  END

( 111 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           T              /            k
           Z           L0>0            R0
           Y             q0            t0
           X              z            D

   ( execution time = 3 seconds )

Example:      L0 = 0.61  ;  q0 =  -0.6   ;   z  = 7

         0.61  ENTER^
        -0.6   ENTER^
           7     XEQ "TOL"  >>>>  D = 1.461832570 1010   light-years
                                      RDN   t0 = 1.556294166 1010   years
                                       RDN  R0 = 2.222433469 1010   light-years
                                       RDN   k  = -1  ( hyperbolic Universe )

-Neglecting matter is probably not realistic, unless an unknown negative density of matter
  exactly compensate for the positive density of matter in the galaxies ( ???? )
 

            f-2)  Radiation+Dust:  Numerical Integration
 
 

Data Registers:   R00 thru R18                          •  R09 = nbre of subintervals if F03 is set  ( unused otherwise )

      R10:  z + 1             R13:   ¶0                                                      R16:   proportional to t0
      R11:  q0                 R14:   2( q0 + L0 - ¶0 )  = (Omega)mat         R17:   proportional to D
      R12:  L0                 R15:  2q0 + 3L0 - ¶0 - 1                              R18:  temp

 -when the program stops:   R00 = k       ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
                                              R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )

Flags:  F00  set flag F00  if you don't want to (re-)calculate the age of the Universe t0
             F03  set flag F03 to fix the number n of subintervals used by "GL3"  and  store n into R09
                              if F03 is clear, n is doubled until 2 consecutive rounded values are equal.
             F10

Subroutine:  "GL3"  Gauss-Legendre 3-point formula  ( cf "Numerical Integration for the HP-41" )
 
 

  01  LBL "COSMO4"
  02  DEG
  03  SF 10
  04  STO 10
  05  "S"
  06  ASTO 00
  07  CLX
  08  STO 01
  09  +
  10  STO 11
  11  X<>Y
  12  STO 12
  13  +
  14  X<>Y
  15  STO 13
  16  -
  17  ST+ X
  18  STO 14
  19  +
  20  RCL 12
  21  +
  22  1
  23  STO 02
  24  ST+ 10
  25  -
  26  STO 15
  27  FC? 00
  28  XEQ 01
  29  FC? 00
  30  STO 16
  31  RCL 10
  32  SQRT
  33  1/X
  34  STO 01
  35  XEQ 01
  36  STO 17
  37  CF 10
  38  XEQ 01
  39  STO 02
  40  STO 04
  41  RCL 15
  42  ABS
  43  SQRT
  44  X=0?
  45  SIGN
  46  1/X
  47  STO 03
  48  STO 06
  49  /
  50  ENTER^
  51  E^X-1
  52  LASTX
  53  CHS
  54  E^X-1
  55  -
  56  2
  57  /
  58  RCL Y
  59  R-D
  60  SIN
  61  RCL 15
  62  X#0?
  63  ENTER^
  64  X>0?
  65  ENTER^
  66  GTO 03
  67  LBL 01
  68  CLX
  69  STO 18
  70  SIGN
  71  FS? 03
  72  RCL 09
  73  STO 03
  74  LBL 02
  75  RCL 03
  76  FC? 03
  77  ST+ 03
  78  XEQ "GL3"
  79  ST+ X
  80  VIEW X
  81  FS? 03
  82  RTN
  83  RND                      If flag F03 is clear, the accuracy is controlled by the display format
  84  ENTER^
  85  X<> 18
  86  X#Y?
  87  GTO 02
  88  LASTX
  89  RTN
  90  LBL "S"
  91  X^2
  92  ENTER^
  93  ENTER^
  94  X^2
  95  RCL 12
  96  *
  97  RCL 15
  98  -
  99  *
100  RCL 14
101  +
102  *
103  RCL 13
104  +
105  /
106  SQRT
107  FS? 10
108  *
109  RTN
110  LBL 03
111  R^
112  RCL 10
113  *
114  ST* 03
115  RCL 15
116  X#0?
117  SIGN
118  STO 00
119  RCL 16
120  STO 05
121  RCL 17
122  STO 01
123  SF 25
124  137 E8
125  ST* 01
126  ST* 02
127  ST* 03
128  ST* 05
129  ST* 06
130  RCL 05
131  SIGN
132  RCL 04
133  RCL 03
134  RCL 02
135  RCL 01
136  CF 25
137  CLA
138  CLD
139  END

( 200 bytes / SIZE 019 )
 
 
      STACK        INPUTS      OUTPUTS
           T             ¶0            VR
           Z             L0            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:     ¶0 = 0.01 ;  L0 = 0.62  ;  q0 =  -0.6   ;   z  = 7

-Let's fix the number of subintervals for the Gaussian integration, for instance:   n = 8   STO 09    SF 03

     0.01  ENTER^
     0.62  ENTER^
     -0.6  ENTER^
        7   XEQ "COSMO4"   >>>>   D  =  1.438258446  1010  l-y  = R01   ( in 80 seconds )
                                            RDN    D0 = 3.419896391  1010  l-y  = R02
                                            RDN    DL = 3.844671484  1011  l-y  = R03
                                            RDN    VR = 2.496274738                 = R04
                                         LASTX     t0 = 1.528106505  1010  y    = R05
                                      [ RCL 06 ]  R0 = 2.315722657  1010  l-y
                                      [ RCL 00 ]   k = -1  ( hyperbolic geometry )

-In this example, almost all digits are significant.
-For our Universe, if we take   ¶0 = 0.000049 ;  L0 = 0.519  ;  q0 =  -0.5   ;   z  = 7 ,  t0 is obtained with only 6 significant digits ( with n = 8 )
-Obviously, the effects of the pressure of radiation are very small and they may be neglected except in the very early age of the Universe.
 

            f-3)  Radiation+Dust:  Elliptic Integrals
 

-There are many subcases and the following program only works if the cosmological constant is positive
  and if the quartic:   L0.y4 + ( 1+¶0-2q0-3L0 ).y2 + 2.( q0+ L0- ¶0 ).y + ¶0
  has 2 distinct non-positive roots and a pair of complex conjugate roots.

-These conditions are precisely satisfied by our Universe!

Formulae:   1°)  The D0 integral is calculated by the following method:

  §ab  [ ( y2+c.y+d ) ( y2+p.y+q ) ] -1/2 dy = 2.RF ( U2 ; U2+T+V ; U2+T-V )

   where  (b-a).U = [ A(a).B(b) ]1/2 + [ A(b).B(a) ]1/2      with   A(y) = ( y2+c.y+d )    and   B(y) = ( y2+p.y+q )
                       T = c.p/2 - d - q
                       V = 2 [ ( c2/4 - d ) ( p2/4 - q ) ]1/2

    2°)  The other integral is of the form   I =  §ab  y [ ( y+m ) ( y+n ) ( y2+p.y+q ) ] -1/2 dy  and involves elliptic integrals of the 3rd kind.

-Unfortunately, I can't have access to the required tables of Elliptic Integrals. So I've split this integral into 2 integrals:

   I =  §ab  ( y+ m ) [ ( y+m ) ( y+n ) ( y2+p.y+q ) ] -1/2 dy  -  §ab  m [ ( y+m ) ( y+n ) ( y2+p.y+q ) ] -1/2 dy

-The integral on the right may be obtained as above.
-And the integral on the left is computed by executing RJZ twice , after the change of variable  u = 1/(y+m)

-But you can certainly simplify the program and reduce execution time if you know the tables of Carlson's elliptic integrals...
 

Data Registers:   R00 thru R34

      R10:  z + 1             R13:   ¶0                                                      R26:   proportional to t0
      R11:  q0                 R14:   2( q0 + L0 - ¶0 )  = (Omega)mat         R29:   proportional to D
      R12:  L0                 R15:  2q0 + 3L0 - ¶0 - 1

 -when the program stops:   R00 = k       ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
                                              R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )

Flags:  F00  set flag F00  if you don't want to re-calculate the age of the Universe t0
              F01 & F02 are used by "P4"

Subroutines:   "P4"  ( cf  "Polynomials for the HP-41" , the new version )
                           "1/Z"  ( cf  "Complex Functions for the HP-41" )
           "RFZ" & "RJZ"   Carlson Elliptic Integrals of the 1st & 3rd kind, complex arguments.
                                     ( cf  "Jacobian elliptic functions and elliptic integrals for the HP-41" )
 

  01  LBL "COSMO5"
  02  STO 10
  03  CLX
  04  +
  05  STO 11
  06  X<>Y
  07  STO 12
  08  +
  09  X<>Y
  10  STO 13
  11  -
  12  ST+ X
  13  STO 14
  14  +
  15  RCL 12
  16  +
  17  1
  18  ST+ 10
  19  -
  20  STO 15
  21  CHS
  22  RCL 14
  23  R^
  24  RCL 12
  25  ST/ T
  26  ST/ Z
  27  /
  28  0
  29  RDN
  30  XEQ "P4"
  31  FS?C 01
  32  FS?C 02
  33  SF 99
  34  STO 33
  35  ST+ X
  36  CHS
  37  STO 18
  38  RDN
  39  STO 34
  40  RDN
  41  X<Y?
  42  X<>Y
  43  STO 16
  44  X<>Y
  45  STO 17
  46  -
  47  STO 32
  48  R^
  49  RCL 33
  50  R-P
  51  X^2
  52  STO 19
  53  RCL 18
  54  +
  55  1
  56  +
  57  STO 21
  58  RCL 10
  59  1/X
  60  RCL 18
  61  -
  62  RCL 10
  63  /
  64  RCL 16
  65  RCL 17
  66  *
  67  STO 20
  68  +
  69  *
  70  SQRT
  71  RCL 10
  72  1/X
  73  RCL 18
  74  +
  75  RCL 10
  76  /
  77  RCL 19
  78  +
  79  RCL 20
  80  RCL 18
  81  -
  82  1
  83  +
  84  STO 22
  85  *
  86  SQRT
  87  +
  88  1
  89  RCL 10
  90  1/X
  91  -
  92  /
  93  X^2
  94  STO 23
  95  RCL 18
  96  X^2
  97  RCL 20
  98  4
  99  *
100  -
101  RCL 19
102  RCL 18
103  2
104  /
105  X^2
106  -
107  *
108  SQRT
109  STO 24
110  RCL 18
111  X^2
112  CHS
113  2
114  /
115  RCL 19
116  -
117  RCL 20
118  -
119  STO 25
120  RCL 23
121  ST+ Y
122  XEQ "RFZ"
123  ST+ X
124  RCL 12
125  SQRT
126  /
127  STO 29
128  FS? 00
129  GTO 01
130  1
131  RCL 34
132  RCL 16
133  ST- Z
134  RCL 33
135  -
136  XEQ "1/Z"
137  STO 30
138  X<>Y
139  STO 31
140  X<>Y
141  RCL 32
142  1/X
143  R^
144  1/X
145  ST+ Y
146  ST+ Z
147  RDN
148  XEQ "RJZ"
149  STO 26
150  STO 27
151  RCL 31
152  RCL 30
153  RCL 32
154  1/X
155  RCL 16
156  CHS
157  X=0?
158  GTO 00
159  1/X
160  ST+ Y
161  ST+ Z
162  RDN
163  XEQ "RJZ"
164  ST- 26
165  LBL 00
166  RCL 33
167  RCL 16
168  -
169  X^2
170  RCL 34
171  X^2
172  +
173  RCL 32
174  *
175  SQRT
176  3
177  *
178  STO 28
179  ST/ 26
180  ST/ 27
181  RCL 24
182  RCL 19
183  RCL 22
184  *
185  SQRT
186  RCL 20
187  RCL 21
188  *
189  SQRT
190  +
191  X^2
192  RCL 25
193  X<>Y
194  ST+ Y
195  XEQ "RFZ"
196  RCL 16
197  *
198  ST+ 26
199  2
200  RCL 12
201  SQRT
202  /
203  ST* 26
204  LBL 01
205  RCL 31
206  RCL 30
207  RCL 10
208  1/X
209  RCL 16
210  -
211  1/X
212  RCL 32
213  1/X
214  X<>Y
215  ST+ Y
216  ST+ Z
217  RDN
218  XEQ "RJZ"
219  RCL 28
220  /
221  RCL 27
222  X<>Y
223  -
224  ST+ X
225  RCL 12
226  SQRT
227  /
228  RCL 29
229  RCL 16
230  *
231  +
232  STO 01
233  RCL 29
234  STO 02
235  STO 04
236  RCL 15
237  ABS
238  SQRT
239  X=0?
240  SIGN
241  1/X
242  STO 03
243  STO 06
244  /
245  ENTER^
246  E^X-1
247  LASTX
248  CHS
249  E^X-1
250  -
251  2
252  /
253  RCL Y
254  R-D
255  SIN
256  RCL 15
257  X#0?
258  ENTER^
259  X>0?
260  ENTER^
261  R^
262  RCL 10
263  *
264  ST* 03
265  RCL 15
266  X#0?
267  SIGN
268  STO 00
269  RCL 26
270  STO 05
271  137 E8
272  ST* 01
273  ST* 02
274  ST* 03
275  ST* 05
276  ST* 06
277  RCL 05
278  SIGN
279  RCL 04
280  RCL 03
281  RCL 02
282  RCL 01
283  END

( 424 bytes / SIZE 035 )
 
 
      STACK        INPUTS      OUTPUTS
           T             ¶0            VR
           Z             L0            DL
           Y             q0            D0
           X              z            D
           L              /            t0

Example:     ¶0 = 0.000049 ;  L0 = 0.519  ;  q0 =  -0.5   ;   z  = 7

     0.000049  ENTER^
       0.519      ENTER^
       -0.5        ENTER^
          7         XEQ "COSMO5"   >>>>   D  =  1.413693452  1010  l-y  = R01
                                                    RDN    D0 = 3.437776327  1010  l-y  = R02
                                                    RDN    DL = 4.219648045  1011  l-y  = R03
                                                    RDN    VR = 2.509325786                 = R04
                                                 LASTX     t0 = 1.565753235  1010  y    = R05
                                              [ RCL 06 ]  R0 = 2.058233710  1010  l-y
                                              [ RCL 00 ]   k = -1  ( hyperbolic geometry )

-Execution time = 187 seconds if flag F00 is clear.
-Execution time =  78  seconds if flag F00 is set.

-The Cosmic Background at  T = 2.736 K leads to  ¶0 = 0.000049  if  1/H0 = 1.37 1010 years
- L0 = 0.519 corresponds to the baryonic matter only.
-With  L0 = 0.519  &  q0 =  -0.5  a value of  ¶0 = 0.001  would still be too great:
  the quartic wouldn't satisfy the required conditions anymore and you'd get  NON EXISTENT  line 33
 

References:      Stamatia Mavridès - "L'Univers relativiste" - Masson  ISBN 2-225-36080-7  ( in French )
                            Jean Heidmann - "Introduction à la cosmologie" - PUF  ( in French )
                            Landau & Lifshitz - "Classical Theory of Fields" - Pergamon Press
                            Feige - "Elliptic Integrals for Cosmological Constant Cosmologies" - Astron. Nachr 313 (1992) 3, 139-163
                            Kantowski, Kao, Thomas - "Distance-Redshift Relations in Inhomogeneous Friedmann-Lemaitre-Robertson-Walker Cosmology" -
                                                                       The Astrophysical Journal, 545: 549-560
 
 
 
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall