hp41programs

General Relativity and Cosmology for the HP-41

General Relativity and Cosmology for the HP-41


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

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

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 numbers are linked by the relation:  ( 1/H0 ) ( years ) = 977.79222143 109 / H0 ( km/s/Mpc )
-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)

        M = Mass ( expressed in Solar Masses )
         T = Orbital Revolution ( in days )
         e  = eccentricity
  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:     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)

       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)

      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  is expressed in Solar 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)

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

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:
-Replace line 08 by 74243 E-32   if   M is expressed in kg  and  r in meter
 
 

 01  LBL "GRZ2"
 02  1/X
 03  CHS
 04  X<>Y
 05  1/X
 06  +
 07  *
 08  2.1208 E-6
 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 loses 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  represents 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 
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 
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 
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: /

-Line 29:  here,  X-register = L0    Add  STO 00  after this line if you want to save this number.
-Line 40:  here,  X-register = a      Add  STO 01  after this line if you want to save this value.
 
 

  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  / 
  30  STO P 
  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 
  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" )

-You'll get  DATA ERROR  line 20  if the density of matter is negative  i-e  q0 + L0 < 0
 
 

  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 
  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 
  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 
  87   E99 
  88  STO 16 
  89  SIGN
  90  STO 03
  91  LBL 03 
  92  RCL 03
  93  ST+ 03 
  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

       FIX 4

     0.7  ENTER^
   -0.6  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.3283 1010 ly         = R01            ( in 93s )             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.

Notes:

-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 )
-Lines 86 to 103  use Gauss Legendre 3-point formula
  the number n of subintervals ( in register R03 ) is doubled until 2 consecutive rounded values are equal
  thus, the accuracy is controlled by the display format

-Actually, for our Universe ( q0 of the order of  -0.5 and L0 of the order of 0.6 ),
  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

-With our Universe, n = 10 produces almost 10 significant figures!

-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 < y  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 
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 
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 
311  1/X 
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

         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" Carlson elliptic integrals of the first and third kind  complex arguments ---------------------------------------------------

-Lines 32 to 69 may be deleted if you don't want to calculate Rmax and  T
-Lines 59 to 64 are only useful if the Universe is an Eddington model ( finite Rmax but no Big Crunch )
 
 

 01  LBL "COSMO3"
  02  1
  03  +
  04  STO 16
  05  FS? 00
  06  GTO 04
  07  CLX
  08  +
  09  + 
  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 
  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? 
  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

       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

       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 yield:

          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:

     t0 = (c/H0)     §01     y. [  (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad  ] -1/2 dy
    D  = (c/H0)  §y(em)1  y. [  (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad  ] -1/2 dy           y(em) = y at the instant of emission
    D0 = (c/H0)  §y(em)1      [  (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad  ] -1/2 dy               z + 1 = R0/R = 1/yem

  where  (Omega)tot = (Omega)mat + (Omega)lambda + (Omega)rad
 

            f-1)  Pure Radiation
 

-Here,  (Omega)mat = 0
-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-¶0-L0)/(2.L0) + 1/L01/2 ] / [ 1/(z+1)2 + (1-¶0-L0)/(2.L0) + ( 1/(z+1)4 + ( (1-¶0-L0)/(z+1)2 + 0 )/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 M
07  SIGN
08  X<>Y
09  STO T
10  -
11  X<>Y
12  STO N
13  -
14  STO O
15  RCL M
16  *
17  RCL N      
18  +
19  X<>Y
20  /
21  RCL M
22  X^2
23  +
24  SQRT
25  RCL M
26  +
27  RCL O
28  R^
29  ST+ X
30  /
31  ST+ Y
32  STO M      
33  R^
34  ST/ N
35  SQRT
36  1/X
37  +
38  1
39  +
40  STO Z
41  X<>Y
42  /
43  LN
44  X<> M
45  RCL N
46  SQRT       
47  +
48  /
49  LN
50  X<>Y
51  SQRT
52  ST+ X
53  ST/ M
54  /
55  RCL O
56  CHS
57  ABS
58  LASTX
59  X#0?
60  SIGN
61  RDN
62  X=0?
63  SIGN
64  SQRT       
65  1/X
66  STO Z 
67  CLX
68  137 E8
69  ST* M
70  ST* Z
71  *
72  RCL M
73  CLA
74  END

 
   ( 109 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           T              /            k
           Z             ¶0            R0
           Y             L0            t0
           X              z            D

   ( execution time = 3 seconds )

Example:       ¶0 =  (Omega)rad = 49 10 -6  ,  L0 = (Omega)lambda = 0.61  ,   z  = 7

      49 E-6  ENTER^
        0.61    ENTER^
           7      XEQ "TOL"  >>>>  D = 1.564459530 1010   light-years
                                        RDN   t0 = 1.814132844 1010   years
                                        RDN  R0 = 2.193893533 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 R17                          •  R09 = nbre of subintervals if F03 is set  ( unused otherwise )

      R10 =  z + 1                  R13 =  (Omega)rad
      R11 =  (Omega)mat       R14  = (Omega)tot - 1
      R12 =  (Omega)lambda  R15-R16-R17:  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" )

-If flag F03 is clear, the accuracy is controlled by the display format
 
 

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

 
   ( 186 bytes / SIZE 018 )
 
 

      STACK        INPUTS      OUTPUTS
           T             ¶0            VR
           Z             L0            DL
           Y      (Omega)mat            D0
           X              z            D
           L             /            t0

 
Example1:      With    (Omega)mat = 0.044 ,  (Omega)lambda  = 0.519 ,  (Omega)rad  = 0.000049   &   z = 7

  •   for instance,  FIX 4    CF 03

     0.000049  ENTER^
        0.519     ENTER^
        0.044     ENTER^
            7        XEQ "COSMO4"  >>>>   D  =  1.4046  E10  l-y     = R01               ---Execution time = 94s---
                                                    RDN    D0 = 3.4012  E10  l-y     = R02
                                                    RDN    DL = 4.1177  E11  l-y     = R03
                                                    RDN    VR = 2.4826                   = R04         ( unit c = 1 )
                                                  LASTX   t0  = 1.5507  E10  years = R05

   and   R00 = k = -1  ( hyperbolic space )
           R06 =  R0 = 2.0725  E10  light-years

-For another redshift - in the same Universe - Set flag F00
 

Example2:      With    (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,  (Omega)rad  = 0.000049  &   z = 7

    ( These Omega-values are close to those suggested by WMAP = NASA's Wilkinson Microwave Anisotropy Probe )

  •   for instance,  SF 03  ,   n = 8  STO 09

     0.000049  ENTER^
        0.732     ENTER^
        0.271     ENTER^
           7         XEQ "COSMO4"  >>>>   D  =  1.2823287  E10  l-y     = R01               ---Execution time = 81s---
                                                    RDN    D0 = 2.8607721  E10  l-y     = R02
                                                    RDN    DL = 2.2835499  E11  l-y     = R03
                                                    RDN    VR = 2.0881548                   = R04         ( unit c = 1 )
                                                  LASTX   t0  = 1.3596702  E10  years = R05

   and   R00 = k = +1  ( spherical space )
           R06 =  R0 = 2.4810862  E11   light-years

-For another redshift - in the same Universe - Set flag F00
 

            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:   (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad

  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 R32

      R10 =  z + 1                  R13 =  (Omega)rad
      R11 =  (Omega)mat       R14  = (Omega)tot - 1
      R12 =  (Omega)lambda  R15 thru R32:  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

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

-Line 124 is a three-byte  GTO 01
 
 

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

 
   ( 416 bytes / SIZE 033 )
 
 

      STACK        INPUTS      OUTPUTS
           T             ¶0            VR
           Z             L0            DL
           Y      (Omega)mat            D0
           X              z            D
           L              /            t0

 
Example1:      With    (Omega)mat = 0.044 ,  (Omega)lambda  = 0.519 ,  (Omega)rad  = 0.000049   &   z = 7

      0.000049  ENTER^
        0.519     ENTER^
        0.044     ENTER^
            7        XEQ "COSMO5"  >>>>   D  =  1.404555850  E10  l-y     = R01               ---Execution time = 63s---
                                                    RDN    D0 = 3.401195936  E10  l-y     = R02
                                                    RDN    DL = 4.117670678  E11  l-y     = R03
                                                    RDN    VR = 2.482624771                   = R04         ( unit c = 1 )
                                                  LASTX   t0  = 1.550654988  E10  years = R05

   and   R00 = k = -1  ( hyperbolic space )
           R06 =  R0 = 2.072546096  E10  light-years

-For another redshift - in the same Universe - Set flag F00
 

Example2:      With    (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,  (Omega)rad  = 0.000049  &   z = 7

    ( These Omega-values are close to those suggested by WMAP = NASA's Wilkinson Microwave Anisotropy Probe )

      0.000049  ENTER^
        0.732     ENTER^
        0.271     ENTER^
           7         XEQ "COSMO5"  >>>>   D  =  1.282328697  E10  l-y     = R01               ---Execution time = 53s---
                                                    RDN    D0 = 2.860772072  E10  l-y     = R02
                                                    RDN    DL = 2.283549897  E11  l-y     = R03
                                                    RDN    VR = 2.088154797                   = R04         ( unit c = 1 )
                                                  LASTX   t0  = 1.359670612  E10  years = R05

   and   R00 = k = +1  ( spherical space )
           R06 =  R0 = 2.481086210  E11   light-years

-For another redshift - in the same Universe - Set flag F00

Notes:

-If F00 is clear, execution time ~  30 seconds.
-These execution times are exact if you have used the M-Code functions  RFZ & RJZ
-Otherwise, they are longer.

-The M-code functions  RFZ & RJZ  also allow to reduce the SIZE, since registers  R00 thru R09 are unused before line 227 after "P4" is executed.

-If the quartic does not satisfy the required conditions, you'll get  NON EXISTENT  line 28
 

References:

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