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  q₀ = 1/2     ( q = deceleration parameter )
     a-2)  Spherical models  q₀ > 1/2
     a-3)  Hyperbolic models  0 < q₀ < 1/2

   b)  Empty Universes:   L₀ + q₀ = 0

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

   c) Eddington-Lemaitre Universes ( critical parameters )

   d) Euclidean Universes

     d-1)   -1 < q₀ < 1/2
     d-2)    q₀ > 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  m³ kg^-1 s^-2
              Heliocentric Gravitational Constant = 1.3271244 10²⁰  m³ s^-2
              Geocentric Gravitational Constant = 3.986004356 10¹⁴  m³ s^-2
              Mass of the Sun = M = 1.98892 10³⁰  kg
              Mass of the Earth = m = 5.97370 10²⁴  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  10¹⁵ m
              1 Megaparsec = 1 Mpc = 3.085677581 10²² m
              Stefan's Constant = a = 7.565656  10 ^-16 J.m³.K ^-1
              Temperature of the Cosmic Background = T = 2.736 K
              Hubble Constant = H₀ ~ 71 km/s/Mpc ~ 2.301 10^-18 s

-In the following programs, I've choosen  1/H₀ = 13.7 10⁹ years which corresponds to  H₀ = 71.37169499 km/s/Mpc = 2.312999111 10^-18 second 

-These numbers are linked by the relation:  ( 1/H₀ ) ( years ) = 977.79222143 10⁹ / H₀ ( 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/m³ , 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:

   ds² = (1-2GM/(c².r)) c².dt² - dr²/(1-2GM/(c².r)) - r²(d(theta)² + sin²(theta) d(phi)² )     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.c²(1-e²)) radians per revolution
        Third Kepler's law:  a³/T² = GM/(4.PI²)

 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:

    m₁ = 1.25 Solar Mass  ;   m₂ = 1.34 Solar Mass   whence   M = m₁ +  m₂ = 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/(c²d)  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:

  d²u/d(phi)² + u = 3.GM/c² u²     where   u = 1/r    and ( r, phi ) = polar coordinates.

-A "funny" solution is  r = Constant = 3.GM/c² = 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/c²/r_O)/(1-2GM/c²/r_S))^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 )
        r_S = distance source-gravitational body     ( in Solar radius )
        r_O = distance observer-gravitational body  ( in Solar radius )

Example:     M = 2 ,  r_S = 3 ,  r_O = 4

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

-If  M is expressed in kg and  r_S ,  r_O are in meters, replace line 02 by  14849 E-31
-When M = mass of the Earth, r_O = radius of the Earth and r_O-r_S = 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/(c²r) are small (  GM/(c²r) << 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 ( m₁ & m₂ ) is slowly decreasing according to the differential equation:

   dr/dt = -64G³ m₁.m₂ ( m₁ + m₂ ) / ( 5.c⁵.r³ )

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

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  m_i are expressed in Solar Masses and  r_j  are expressed in Astronomical Units
  >>>   The result T is expressed in years

Example:   Let's take once again the binary pulsar PSR J0737-3039      m₁ = 1.25  ,   m₂ = 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 m_i are in kg and  r_j are in meters, replace line 14 with  1.592026 E71
 

3°) Cosmology
 

-Einstein's equations   S_ij - (1/2) g_ij ( S + 2.Lambda ) = ( 8.PI.G/c⁴ ) T_ij   are now applied to a homogeneous and isotropic Universe.
    (  S_ij is the tensor of curvature , S = S_ijS^ij , Lambda is the cosmological constant and  T_ij describes the content of the Universe  )

-This "cosmological principle" leads to the Robertson-Walker metric:

    ds² = c²dt² - R²(t)/(1-kr²/4)² . ( dr² + r²d(theta)² + r²sin²(theta)d(phi)² )

      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,  T_ij = ( (rho).c² + p ) u_iu_j - p g_ij      ( g_ij is the metric tensor ; rho = density of matter and energy ; p = pressure ; u_i = dx_i/ds = 4-velocity )
  and Einstein's equations yield:

  2 R(t).d²R/dt² + (dR/dt)² + k.c² = ( -(8.PI.G/c² ).p + (Lambda).c² ).R²(t)
     (dR/dt)² + k.c²  = ( (8.PI.G/3) (rho) + (Lambda/3).c² ).R²(t)

_________________________________________________________________________________________________________________________

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

   With  rho = 5 10 ^-28 kg/m³ it yields   R = 4.63 10²⁶ m ~ 15 Gpc  ,  Lambda = 4.665 10 ^-54 m^-2
   Einstein's Universe is a 3-sphere and its volume is  2.PI².R³ = 1.959 10⁸¹ m³ ~ 66700 Gpc³
   Its mass is M = 9.796 10⁵³ 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.(d²R/dt²)/(dR/dt)²   =  deceleration parameter
  (Omega)_mat = 8.PI.G.(rho)/(3H²) =  density parameter
  (Omega)_vac = (Lambda.c²)/(3H²) =  L

-The equations reduce to:

          2.q + 2.L = 8.PI.G.(rho)/(3H²)
     2.q + 3.L - 1 = k.c²/(H²R²)

  and lead to the integral      H₀.t = §  [ ( 2q₀+2.L₀ )/y + ( 1-2q₀-3L₀ ) + L₀.y² ] ^-1/2 dy

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

-Knowing the cosmological redshift z of a galaxy, and the values of the current parameters q₀ , L₀
 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 )
             D₀ = comoving distance ( in light-years )
             D_L = luminosity-distance ( in light-years )
             V_R = recessional velocity ( speed of light = 1 )
              t₀  = age of the Universe ( in years )

     D  = (c/H₀)  §_y(em)¹  [ ( 2q₀+2.L₀ )/y + ( 1-2q₀-3L₀ ) + L₀.y² ] ^-1/2 dy               y(em) = y at the instant of emission     z + 1 = R₀/R = 1/y_em
     D₀ = (c/H₀)  §_y(em)¹  [ ( 2q₀+2.L₀ ).y + ( 1-2q₀-3L₀ ).y² + L₀.y⁴ ] ^-1/2 dy           D₀ = D_now

                                                                F(x) = Sinh(x)  if  k = -1
     D_L = R₀ ( z + 1 ) F(D₀/R₀)    where    F(x) =   x         if  k = 0
                                                               F(x) = Sin(x)    if  k = +1

     V_R = H₀.D₀
      t₀  =  (1/H₀)  §₀¹  [ ( 2q₀+2.L₀ )/y + ( 1-2q₀-3L₀ ) + L₀.y² ] ^-1/2 dy

-We also have  k = sign( 2q₀+3L₀-1 )   and   R₀ = (c/H₀) k.( 2q₀+3L₀-1 ) ^-1/2   R₀ may be choosen arbitrarily ( positive ) if k = 0

-Since  z + 1 = R₀/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                C₁
     *                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
                                            |                                       *                          C₂
                                                                                        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  & EAC₁  there is a Big Bang and a Big Crunch

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

-Domains3&4:  OAB  & C₁ABC₂  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:    WBC₂  there is no Big Bang and no Big Crunch. R(t) has a positive minimum ( bounce models )

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

-Critical parameters:

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

       R(t)
          |                                                                  *
           |                                                         *                 Einstein-Lemaitre models:   Line )BC₂)      There is no Big Bang and  R(t)  = R_E  for t = - infinity
           |                                           *
           |                          *
           |*
           |---------------------------------------------------------------------------- Einstein's model   R(t) = R_E = constant
           |                                                                       *
           |                           *
           |           *
           |    *                                                                      Einstein-Eddington models:  Line )AC₁)      There is a Big Bang and R(t)  = R_E  for t = + infinity
     ------|*-------------------------------------------------------------------------- t
 
 

_________________________________________________________________________________________________________________________
 
 

        a) Friedmann Universes: No Cosmological Constant , L₀ = 0

         a-1)  Einstein-de Sitter Model:  q₀ = 1/2 ; L₀ = 0

-This is an Euclidean Universe:   2q₀ + 3L₀ - 1 = 0     whence  k = 0

Formulae:

    D = (c/H₀).(2/3).( 1 - (1 + z )^-3/2 )
   D₀ = 2(c/H₀).( 1 - (1 + z )^-1/2 )
   D_L = D₀ ( 1 + z )                               R(t) = (3H₀t/2)^2/3    if we choose  R₀ = 1
   V_R = H₀ D₀
    t₀ =  (2/3)H₀
 
 

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 10⁹ ly
                            RDN   D₀ = 17.713 10⁹ ly
                            RDN   D_L = 141.701 10⁹ ly
                            RDN   V_R =  1.2929
                         LASTX  t₀  =  9.133 10⁹ years
 

         a-2)  Spherical Models:  q₀ > 1/2 ; L₀ = 0
 

-Here,  k = +1

Formulae:

    D = (c/H₀/(2q₀-1)).[ -1 + (2q₀z+1)^1/2/(z+1) - (q₀/(2q₀-1)^1/2) ( Arccos ((2q₀-1)/q₀-1) - Arccos ((2q₀-1)/(q₀(z+1)) -1) ]
   D₀ = (c/H₀/(2q₀-1)^1/2).[ - Arccos ((2q₀-1)/q₀-1) + Arccos ((2q₀-1)/(q₀(z+1)) -1) ]
   D_L = R₀ ( 1 + z ) Sin D₀/R₀                                                                                                R₀ = c/H₀/(2q₀-1)^1/2
   V_R = H₀ D₀
    t₀ =  (1/H₀/(2q₀-1)) [ -1 - (q₀/(2q₀-1)^1/2) ( Arccos ((2q₀-1)/q₀-1) - PI ]

-The period of these Universes ( time between a Big Bang and a Big Crunch ) is  T = (1/H₀)(2.pi.q₀).(2q₀-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:      q₀ = 2  ;  z = 3

    2  ENTER^
    3  XEQ "FRSPH"  >>>>  D  =  5.871 10⁹ ly
                                 RDN   D₀ =  9.482 10⁹ ly
                                 RDN   D_L = 29.474 10⁹ ly
                                 RDN   V_R =  0.6921
                               LASTX  t₀  =  6.477 10⁹ years

Note:

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

         a-3)  Hyperbolic Models:  0 < q₀ < 1/2 ; L₀ = 0
 

-Here,  k = -1

Formulae:

    D = (c/H₀/(1-2q₀)).[ 1 - (2q₀z+1)^1/2/(z+1) - (q₀/(1-2q₀)^1/2) ( Arccosh ((1-q₀)/q₀) - Arccos ((1-2q₀)/(q₀(z+1)) +1) ]
   D₀ = (c/H₀/(1-2q₀)^1/2).[  Arccosh ((1-q₀)/q₀) - Arccosh ((1-2q₀)/(q₀(z+1)) +1) ]
   D_L = R₀ ( 1 + z ) Sinh D₀/R₀                                                                                          R₀ = c/H₀/(1-2q₀)^1/2
   V_R = H₀ D₀
    t₀ =  (1/H₀/(1-2q₀)) [ 1 - (q₀/(1-2q₀)^1/2)  Arccosh ((1-q₀)/q₀)  ]
 
 

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:      q₀ = 0.4  ;  z = 7

   0.4  ENTER^
    7  XEQ "FRHYP"  >>>>  D  =   9.087 10⁹ ly
                                 RDN   D₀ = 18.708 10⁹ ly
                                 RDN   D_L = 159.140 10⁹ ly
                                 RDN   V_R =  1.3655
                               LASTX  t₀  =  9.534 10⁹ years
 

        b) Empty Universes:   q₀ + L₀ = 0
 

         b-1)  Milne Model:   q₀ = 0
 

  k = -1  ( hyperbolic geometry )

Formulae:

    D = (c/H₀).z/(z+1)
   D₀ = (c/H₀).Ln(z+1)
   D_L = R₀ z(z+2)/2           R₀ = c/H₀        R(t) = c.t
   V_R = H₀ D₀
    t₀ =  1/H₀
 
 

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 10⁹ ly
                              RDN   D₀ = 28.488 10⁹ ly
                              RDN   D_L = 431.550 10⁹ ly
                              RDN   V_R =  2.0794
                             LASTX  t₀  = 13.700 10⁹ years

         b-2) de Sitter  Model:   q₀ = -1
 

  k = 0  ( Euclidean geometry )

Formulae:

    D = (c/H₀).Ln(z+1)
   D₀ = (c/H₀).z
   D_L = (c/H₀).z(z+1)
   V_R = H₀ D₀
    t₀ =  infinity              R(t) = exp(H₀.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 10⁹ ly
                           RDN   D₀ =  95.900 10⁹ ly
                           RDN   D_L = 767.200 10⁹ ly
                           RDN   V_R =  7.0000
                         LASTX  t₀  =  9.9999 10⁹⁹ years = infinity
 

         b-3)  -1  <  q₀  <  0
 

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

Formulae:       k = -1   ( hyperbolic geometry )

    D = (c/H₀/(-q₀)^1/2). [ Arcsinh (1/b) - Arcsinh (1/b/(z+1)) ]      where   b = ((1+q₀)/(-q₀))^1/2
   D₀ = (c/H₀/(1+q₀)^1/2). Ln ( ( 1/b² + (z+1)² )^1/2 + z+1 )/((1/b²+1)^1/2+1 )
   D_L = R₀ ( 1 + z ) Sinh D₀/R₀          R₀ = c/H₀/(1+q₀)^1/2       R(t) = R₀.b Sinh ( H₀.t.(-q₀)^1/2 )
   V_R = H₀ D₀
    t₀ =  1/H₀/(-q₀)^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:      q₀ = -0.3 ;  z = 7

 -0.3  ENTER^
    7   XEQ "VAC"   >>>>  D  =  13.341 10⁹ ly
                                RDN   D₀ =  32.552 10⁹ ly
                                RDN   D_L = 469.215 10⁹ ly
                                RDN   V_R =  2.3761
                              LASTX  t₀  =  15.386 10⁹ years 

-We can save several bytes if we don't use status registers M, N, O ....:

01 LBL "VAC"  02 1  03 STO 00  04 STO 06  05 +  06 STO 01  07 STO 04            08 X<>Y  09 ENTER  10 ST+ 06  11 CHS  12 STO 02  13 LASTX  14 -

15 /  16 STO 03  17 SQRT  18 STO 05  19 LASTX  20 RCL 00            21 +  22 SQRT  23 ST+ 00  24 +  25 ENTER  26 R^  27 ST* Y  28 X^2

29 RCL 03  30 +  31 SQRT  32 ST+ 04  33 RCL 05            34 +  35 /  36 LN  37 X<> 04  38 RCL 00  39 /  40 LN  41 LASTX  42 ENTER

43 1/X  44 -  45 RCL 01            46 *  47 2  48 /  49 R^  50 LN  51 RCL 02  52 SQRT  53 ST/ 04  54 /  55 RCL 06  56 SQRT

57 ST/ Z  58 ST/ T  59 RDN  60 SIGN  61 CLX  62 RCL Z  63 137 E8  64 ST* L  65 ST* 04  66 ST* Y  67 ST* Z  68 X<> 04            69 END

 
    ( 95 bytes / SIZE 007 )
 
 

STACK	INPUTS	OUTPUTS
T	/	VR
Z	/	DL
Y	q0	D0
X	z	D
L	/	t0

 
Example:      q₀ = -0.3 ;  z = 7

 -0.3  ENTER^
    7   XEQ "VAC"   >>>>  D  =  13.341 10⁹ ly
                                RDN   D₀ =  32.552 10⁹ ly
                                RDN   D_L = 469.215 10⁹ ly
                                RDN   V_R =  2.3761
                              LASTX  t₀  =  15.386 10⁹ years 

-With  13.77172143 ( i-e H₀ = 71 km/s/Mpc )  instead of  137 E8 ( line 63 ), the results are:

 -0.3  ENTER^
    7   XEQ "VAC"   >>>>  D  =  13.41113179 Gly
                                RDN   D₀ =  32.72273316 Gly
                                RDN   D_L = 471.6710452 Gly
                                RDN   V_R =  2.376081547
                              LASTX  t₀  =  15.46638674 Gy
 

         b-4)  q₀  >  0
 

   k = -1   ( hyperbolic geometry )

Formulae:

    D = (c/H₀/q₀^1/2). [ Arcsin (1/b) - Arcsin (1/b/(z+1)) ]      where   b = ((1+q₀)/q₀)^1/2
   D₀ = (c/H₀/(1+q₀)^1/2). Ln ( ( ( -1/b² + (z+1)² )^1/2 + z+1 )/((-1/b²+1)^1/2+1 )
   D_L = R₀ ( 1 + z ) Sinh D₀/R₀          R₀ = c/H₀/(1+q₀)^1/2       R(t) = R₀.b Sin ( H₀.t.q₀^1/2 )
   V_R = H₀ D₀
    t₀ =  1/H₀/q₀^1/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:      q₀ =  0.3 ;  z = 7

   0.3  ENTER^
    7   XEQ "VAC2"   >>>>  D  =  11.031 10⁹ ly
                                 RDN   D₀ =  25.737 10⁹ ly
                                 RDN   D_L = 403.673 10⁹ ly
                                 RDN   V_R =  1.8786
                               LASTX  t₀  =  12.534 10⁹ years
 

         b-5)  q₀  <  -1
 

    k = +1  ( spherical geometry )

Formulae:

    D = (c/H₀/(-q₀)^1/2). [ Arccosh (1/b) - Arccosh (1/b/(z+1)) ]      where   b = ((1+q₀)/q₀)^1/2
   D₀ = (c/H₀/(-1-q₀)^1/2). [ Arccos b - Arccos b(z+1) ]
   D_L = R₀ ( 1 + z ) Sin D₀/R₀          R₀ = c/H₀/(-1-q₀)^1/2       R(t) = R₀.b Cosh ( H₀.t.(-q₀)^1/2 )      R_min =  R₀.b
   V_R = H₀ D₀
    t₀ =  1/H₀/(-q₀)^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:      q₀ =  -1.1 ;  z = 2

  -1.1  ENTER^
     2   XEQ "VAC3"   >>>>  D  =  18.458 10⁹ ly
                                  RDN   D₀ =  35.699 10⁹ ly
                                 RDN   D_L =  95.381 10⁹ ly
                                  RDN   V_R =  2.6057
                                LASTX  t₀  =  24.408 10⁹ years

***  In these models, there is no Big Bang:

  t₀ = 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:   z_max = (q₀/(q₀+1))^1/2 - 1
-With  q₀ = -1.1  ,  z_max = 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,  ( 3L₀+2q₀-1)³ = 27 L₀(q₀+L₀)²   and  L₀ is found by:    6.L₀ = ( q₀ + 1 )² - 6.( q₀ + 1 ) + 6 + [ ( q₀ + 1 )⁴ -(4/3).( q₀ + 1 )³ ]^1/2

  Formulae:         with  a = (1/2).(1+q₀/L₀) ^-1/3

  H₀.D/c = (1/(2a)) sign(-q₀).[2a(q₀+L₀)] ^-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-a^1/2 ).( (1+a)^1/2+a^1/2 )²/( (z+1+a)^1/2+a^1/2 ) ] }

 H₀.D₀/c =  sign(-q₀).[2a(q₀+L₀)] ^-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 ) ]

  H₀.t₀ = (1/(2a)) [2a(q₀+L₀)] ^-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+a^1/2 )² }  if  q₀ > 1/2
       t₀ = infinity  if  q₀ < -1

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

Data Registers: /
Flag:  F10
Subroutines: /

-Line 29:  here,  X-register = L₀    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  q₀ > 1/2   ( Eddington models )
     or   q₀ < -1     ( Lemaitre models )

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

Example1:      q₀ =  -1.1 ;  z = 2      whence   L₀ = 1.107976565  ( line 29 )

  -1.1  ENTER^
     2   XEQ "LEM"    >>>>  D  =  17.355 10⁹ ly
                                  RDN   D₀ =  32.783 10⁹ ly
                                  RDN   D_L =  87.124 10⁹ ly
                                  RDN   V_R =  2.3930
                                LASTX  t₀  =  9.9999 10⁹⁹ years = infinity

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

Example2:      q₀ =  2.375 ;  z = 2      whence   L₀ = 1  ( line 29 )

   2.375  ENTER^
     2   XEQ "LEM"    >>>>  D  =   5.0238 10⁹ ly
                                  RDN   D₀ =  7.4017 10⁹ ly
                                 RDN   D_L =  15.599 10⁹ ly
                                  RDN   V_R =  0.5403
                                LASTX  t₀  =  5.7825 10⁹ years

-If  q₀ > 1/2  ,  y_max = 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.L₀ + 2.q₀ - 1 = 0

   q₀ = 1/2  is Einstein-de Sitter model
   q₀ = -1   is de Sitter model

-For  other q₀-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 < q₀ < 1/2
 

Formulae:    H₀.t₀ =  (3-6q₀) ^-1/2  Ln ( (1+a)^1/2 + a^1/2 )²         with  a = (1-2q₀)/(2+2q₀)

    H₀.D/c = (3-6q₀) ^-1/2  Ln [ ( ( (z+1)³ + a )^1/2 - a^1/2 ) ( ( (z+1)³ + a )^1/2 + a^1/2 ) ^-1 ( (1+a)^1/2 + a^1/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:     q₀ = -0.7  &  z = 7      ( whence  L₀ = 0.8 )

    -0.7  ENTER^
       7   XEQ "EUCL"   >>>>   D = 1.3840  10¹⁰ ly
                                    X<>Y   t₀ = 1.4742  10¹⁰ years
 

           d-2)  q₀ > 1/2
 

Formulae:    H₀.t₀ =  (-3+6q₀) ^-1/2 [ PI - 2.Arctan ( (1+a)^1/2  (-a) ^-1/2 ) ]      with  a = (1-2q₀)/(2+2q₀)

    H₀.D/c = 2.(-3+6q₀) ^-1/2   [ Arctan ( ( (z+1)³ + 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:     q₀ = 2  &  z = 7      ( whence  L₀ = -1 )

       2   ENTER^
       7   XEQ "EUCL"   >>>>   D = 6.8878  10⁹ ly
                                    X<>Y   t₀ = 7.1733  10⁹ years
 

        e) More General Models
 

           e-1)  Numerical Integration ( program#1 )
 

-Three integrals are needed to compute  t₀ , D and D₀
-The following program uses the change of variable  y = u² and Gauss-Legendre 3-point formula.
 

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

             R09:  z + 1               R12:   2(q₀ + L₀)                 R15-R16: temp
             R10:  q₀                   R13:   2q₀ + 3L₀ - 1
             R11:  L₀                   R14:   (1/H₀).t₀

 when the program stops:

    R00 = k     ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
    R01 = D ,  R02 = D₀ , R03 = D_L , R04 = V_R , R05 = t₀ , R06 = R₀ ( current scale factor )

Flags:  F00  set flag F00  if you don't want to re-calculate the age of the Universe t₀
            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  q₀ + L₀ < 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:     L₀ = 0.7  ;  q₀ =  -0.6   ;   z = 7

*** These values correspond to  (Omega)_vac = L₀ = 0.7  &  (Omega)_mat =  2(L₀+q₀) = 0.2   whence (rho)_mat = 3H₀²(L₀+q₀)/(4PI.G) ~ 1.91 10^-27 kg/m³

       FIX 4

     0.7  ENTER^
   -0.6  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.3283 10¹⁰ ly         = R01            ( in 93s )             The successive approximations
                                      RDN   D₀ =  3.0246 10¹⁰ ly         = R02                                        of the 3 integrals (  for  t₀  D₀  D_L )
                                      RDN   D_L =  2.6210 10¹¹ ly         = R03                                        are displayed during the calculations.
                                      RDN   V_R =  2.2077                     = R04
                                    LASTX  t₀  =  1.4168 10¹⁰ years    = R05

-We also have    R06 = 4.3323 10¹⁰ = current scale factor = R(t₀) = R₀
            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)³ >= 27 L(q+L)²   and  q <= -1  there is no Big Bang  and  t₀ = 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 ( q₀ of the order of  -0.5 and L₀ 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 H₀.( t₀ - t )  only )

 for example:    z = 7  means  y = 1/8  (  z = 1/y - 1 ) and  light-travel time = 1.3283 10¹⁰ =  t₀ - t
  whence  t = 8.85 10¹⁰  or  H₀.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:    q₀ =  -0.6  ;  L₀ = 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)/R₀
        |
    3  |                                                                   *
        |
        |                                                               *
        |
    2  |                                                        *
        |                                                    *
        |                                               *
        |                                         *
    1  |                                  * now  ( y = 1 )
        |                         *
        |                * Infl
        |     *
        |*----------|----------|----------|----------|-------  H₀.t
       0              0.5             1            1.5            2

 Here,  y_I = ( 1 + q₀/L₀)^1/3 = 0.522757959      I  is a point of inflexion:     q > 0  for  y < y_I   and  q < 0  for y > y_I
 

Example2:     L₀ = 2/3  ;  q₀ =  -0.5   ;   z  = 7

*** These values correspond to  (Omega)_vac = L₀ = 2/3  &  (Omega)_mat =  2(L₀+q₀) = 1/3   whence (rho)_mat = 3H₀²(L₀+q₀)/(4PI.G) ~ 3.19 10^-27 kg/m³

      2    ENTER^  3  /
   -0.5  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.2123 10¹⁰ ly         = R01            ( in 97s in FIX 4 )
                                      RDN   D₀ =  2.6605 10¹⁰ ly         = R02
                                      RDN   D_L =  2.1284 10¹¹ ly         = R03
                                      RDN   V_R =  1.9420                     = R04
                                    LASTX  t₀  =  1.2822 10¹⁰ years    = R05

-We also have    R06 = 1.3700 10¹⁰ = current scale factor = R(t₀) = R₀
            and        R00 = k = 0 , we have an Euclidean Universe.
 

Example3:     L₀ = 0.6  ;  q₀ =  -0.3   ;   z  = 7

*** These values correspond to  (Omega)_vac = L₀ = 0.4  &  (Omega)_mat =  2(L₀+q₀) = 0.6   whence (rho)_mat = 3H₀²(L₀+q₀)/(4PI.G) ~ 5.74 10^-27 kg/m³

     0.6  ENTER^
   -0.3  ENTER^
     7    XEQ "COSMO"  >>>>  D  =  1.0689 10¹⁰ ly         = R01            ( in 120s in FIX 4 )
                                      RDN   D₀ =  2.2467 10¹⁰ ly         = R02
                                      RDN   D_L =  1.6405 10¹¹ ly         = R03
                                      RDN   V_R =  1.6399                     = R04
                                    LASTX  t₀  =  1.1217 10¹⁰ years    = R05

-We also have    R06 = 3.0634 10¹⁰ = current scale factor = R(t₀) = R₀
            and        R00 = k = +1 , the Universe is spherical.
 

Notes:

-If there is no Big Bang,  the minimum  y_min = R_min/R₀ may be found by solving the cubic equation:   L₀.y³ + ( 1-2q₀-3L₀ ).y +  ( 2q₀+2.L₀ )  = 0

  For example, if    q₀ = -1.1 & L₀ = 1.102   the cubic  1.102 y³ - 0.106 y + 0.004 = 0  has 3 real roots:  0.289201975 ; 0.038320885 ; -0.327522859
  and   y_min = 0.289201975  which corresponds to  z_max = 1/y_min - 1 = 2.457791053

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

  1.102  ENTER^
   -1.1   ENTER^
    2.45  XEQ "COSMO"  yields  D =  2.459 10¹⁰                 in  1h16mn      for only 2 integrals !
                                                   D₀ = 5.599 10¹⁰
                                                   D_L = 1.410 10¹¹
                                                   V_R = 4.087           R₀ = 4.208 10¹⁰

-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 = y_min cosh² v
   b) If  y  has a finite maximum ( oscillating models )  we set:   y = y_max sin² 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(q₀ + L₀)                   R16:   y_max                               R19:  proportional to D
      R11:  q₀                 R14:   -2q₀ - 3L₀ + 1            R17:   proportional to t₀            R20-R21-R22: temp
      R12:  L₀                 R15:  y_min                             R18:   proportional to T

 when the program stops:

      R00 = k       ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
      R01 = D ,  R02 = D₀ , R03 = D_L , R04 = V_R , R05 = t₀ , R06 = R₀ ( current scale factor )
      R07 = R_min  ,  R08 = R_max  ,   R09 = T = Period of an oscillating Universe ( 0 otherwise )

-Here, for bounce Universes,  t₀ 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 t₀and 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:     L₀ = 1.102  ;  q₀ =  -1.1   ;   z  = 2.45

    1.102    ENTER^
    -1.1      ENTER^
     2.45    XEQ "COSMO2"  >>>>  D  =  2.4589 10¹⁰ ly           = R01            ( in 110s in FIX 3 instead of more than 1 hour with "COSMO" )
                                            RDN   D₀ =  5.5990 10¹⁰ ly           = R02
                                             RDN   D_L =  1.4101 10¹¹ ly           = R03
                                             RDN   V_R =  4.0869                       = R04
                                             LASTX  t₀  =  2.5501 10¹⁰ years    = R05

-We also have

   R00 = k = 1   ( Spherical Universe )
   R06 = 4.2079 10¹⁰ = current scale factor = R(t₀) = R₀
   R07 = 1.2169 10¹⁰ = R_min  ( bounce model )
   R08 = 9.9999 10⁹⁹ = R_max = infinity
   R09 = 0  ( non-oscillating Universe )
 

Example2:     L₀ = 0.2  ;  q₀ =  2   ;   z  = 7

     0.2  ENTER^
       2   ENTER^
       7     R/S  >>>>  D  =  6.1763 10⁹ ly          = R01            ( in 158s in FIX 3 )
                     RDN   D₀ =  1.1392 10¹⁰ ly         = R02
                     RDN   D_L =  5.7763 10¹⁰ ly         = R03
                     RDN   V_R =  0.8315                     = R04
                   LASTX  t₀  =  6.3750 10⁹ years     = R05

-We also have

   R00 = k = +1  ( Spherical Universe )
   R06 = 7.2205 10⁹ = current scale factor = R(t₀) = R₀
   R07 = 0 = R_min
   R08 = 9.8405 10⁹ = R_max
   R09 = 3.5658 10¹⁰ = T = time between a Big Bang and a Big Crunch.
 

Example3:     L₀ = 1.108  ;  q₀ =  -1.1   ;   z  = 7

    1.108    ENTER^
    -1.1      ENTER^
       7          R/S   >>>>  D  =  6.7540 10¹⁰ ly           = R01            ( in 11mn25s in FIX 3  )
                          RDN   D₀ =  2.8508 10¹¹ ly           = R02
                          RDN   D_L =  2.6910 10¹¹ ly           = R03
                          RDN   V_R =  20.8088                     = R04
                          LASTX  t₀  =  7.2650 10¹⁰ years    = R05

         R00 = k = 1   ( Spherical Universe )
         R06 = 3.8905 10¹⁰ = current scale factor = R(t₀) = R₀
         R07 = 0  = R_min
         R08 = 9.9999 10⁹⁹ = R_max = infinity
         R09 = 0  ( non-oscillating Universe )

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

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

-  t₀ , D and D₀ 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/H₀)    §₁^z+1 du / u / ( 2(q₀+L₀ ).u³ + (1-2q₀-3L₀ ).u² + L₀ )^1/2         u = 1/y      ;     t₀ corresponds to z = infinity
              D₀ = (c/H₀)    §₁^z+1 du / ( 2(q₀+L₀ ).u³ + (1-2q₀-3L₀ ).u² + L₀ )^1/2

-If the cubic    2(q₀+L₀ ).u³ + (1-2q₀-3L₀ ).u² + L₀ = 0   has 3 real roots  a , b , c  we have:

       D = (c/3/H₀) (2/(q₀ + L₀ ))^1/2  [  R _J ( -a+1 , -b+1 , -c+1 , 1 )  -  R _J ( -a+z+1 , -b+z+1 , -c+z+1 , z+1 )  ]
       D₀ = (c/H₀) (2/(q₀ + L₀ ))^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(q₀+L₀ ).u³ + (1-2q₀-3L₀ ).u² + L₀ = 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/H₀) (2/(q₀ + L₀ ))^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 )  ]
      D₀ = (c/H₀) (2/(q₀ + L₀ ))^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 t₀
             R14:  -b                   R17:   q₀ + L₀                     R20:  y_max
             R15:  -c                   R18:   2q₀ + 3L₀ - 1            R21:  proportional to T

 when the program stops:

     R00 = k         ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
     R01 = D ,  R02 = D₀ , R03 = D_L , R04 = V_R , R05 = t₀ , R06 = R₀ ( current scale factor )
     R07 = R_max  ,  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 t₀and 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 R_max and  T
-Lines 59 to 64 are only useful if the Universe is an Eddington model ( finite R_max 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:     L₀ = 0.7  ;  q₀ =  -0.6   ;   z  = 7

     0.7  ENTER^
   -0.6  ENTER^
      7   XEQ "COSMO3"  >>>>  D  =  1.328255073 10¹⁰ ly         = R01            ( in 140s )      flag F01 is set
                                        RDN   D₀ =  3.024557313 10¹⁰ ly         = R02
                                        RDN   D_L =  2.621046346 10¹¹ ly         = R03
                                        RDN   V_R =  2.207706068                     = R04
                                      LASTX  t₀  =  1.416778742 10¹⁰ years    = R05

       R00 = k = -1 , meaning the Universe is hyperbolic.
       R06 = 4.332320394 10¹⁰ = current scale factor = R(t₀) = R₀
       R07 = 9.999999999 10⁹⁹ = infinity = R_max
       R08 = 0  ( non oscillating Universe )
 

Example2:     L₀ = 0.2  ;  q₀ =  2   ;   z  = 7

     0.2  ENTER^
       2   ENTER^
       7   XEQ "COSMO3"  >>>>  D  =  6.176349709 10⁹ ly          = R01            ( in 145s )      flag F01 is clear
                                        RDN   D₀ =  1.139222269 10¹⁰ ly         = R02
                                        RDN   D_L =  5.776287398 10¹⁰ ly         = R03
                                        RDN   V_R =  0.831549101                     = R04
                                      LASTX  t₀  =  6.375024998 10⁹ years     = R05

       R00 = k = +1  ( Spherical Universe )
       R06 = 7.220533991 10⁹ = current scale factor = R(t₀) = R₀
       R07 = 9.840506634 10⁹ = R_max
       R08 = 3.565763808 10¹⁰ = 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  L₀ = 1.108  &  q₀ =  -1.1  "COSMO3"  yields  t₀  = 7.264885687 10¹⁰ 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)/(3H²) + ¶                 where      ¶ = 8.PI.G.p_rad/(c²H²)  is the pressure parameter = (Omega)_rad
     2.q + 3.L - 1 = k.c²/(H²R²) + ¶                             and        p_rad is the pressure of radiation

   (rho) = (rho)_mat + 3.p_rad/c²   and according to Stefan's law:   3p_rad = (rho)_rad c² = a.T⁴
                                                                                   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 R³ = constant  and   p_rad R⁴ = constant.
-The integrals become:

     t₀ = (1/H₀)     §₀¹     y. [  (Omega)_lambda.y⁴ + ( 1-(Omega)_tot ).y² + (Omega)_mat.y + (Omega)_rad  ] ^-1/2 dy
    D  = (c/H₀)  §_y(em)¹  y. [  (Omega)_lambda.y⁴ + ( 1-(Omega)_tot ).y² + (Omega)_mat.y + (Omega)_rad  ] ^-1/2 dy           y(em) = y at the instant of emission
    D₀ = (c/H₀)  §_y(em)¹      [  (Omega)_lambda.y⁴ + ( 1-(Omega)_tot ).y² + (Omega)_mat.y + (Omega)_rad  ] ^-1/2 dy               z + 1 = R₀/R = 1/y_em

  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 t₀ , the current radius of the Universe R₀ and the constant  k
  ( D & t₀ may be expressed in terms of elementary functions )

Formula:

  H₀.(t₀ - t_e ) = 1/(2.L₀^1/2) Ln { [ 1 + (1-¶₀-L₀)/(2.L₀) + 1/L₀^1/2 ] / [ 1/(z+1)² + (1-¶₀-L₀)/(2.L₀) + ( 1/(z+1)⁴ + ( (1-¶₀-L₀)/(z+1)² + ¶₀ )/L₀ )^1/2 ] }

  ( replace z by +infinity to get  t₀ )

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:       ¶₀ =  (Omega)_rad = 49 10 ^-6  ,  L₀ = (Omega)_lambda = 0.61  ,   z  = 7

      49 E-6  ENTER^
        0.61    ENTER^
           7      XEQ "TOL"  >>>>  D = 1.564459530 10¹⁰   light-years
                                        RDN   t₀ = 1.814132844 10¹⁰   years
                                        RDN  R₀ = 2.193893533 10¹⁰   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 = D₀ , R03 = D_L , R04 = V_R , R05 = t₀ , R06 = R₀ ( current scale factor )

Flags:  F00  set flag F00  if you don't want to (re-)calculate the age of the Universe t₀

             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    D₀ = 3.4012  E10  l-y     = R02
                                                    RDN    D_L = 4.1177  E11  l-y     = R03
                                                    RDN    V_R = 2.4826                   = R04         ( unit c = 1 )
                                                  LASTX   t₀  = 1.5507  E10  years = R05

   and   R00 = k = -1  ( hyperbolic space )
           R06 =  R₀ = 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    D₀ = 2.8607721  E10  l-y     = R02
                                                    RDN    D_L = 2.2835499  E11  l-y     = R03
                                                    RDN    V_R = 2.0881548                   = R04         ( unit c = 1 )
                                                  LASTX   t₀  = 1.3596702  E10  years = R05

   and   R00 = k = +1  ( spherical space )
           R06 =  R₀ = 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.y⁴ + ( 1-(Omega)_tot ).y² + (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 D₀-integral is calculated by the following method:

           §_a^b  [ ( y²+c.y+d ) ( y²+p.y+q ) ] ^-1/2 dy = 2.R_F ( U² ; U²+T+V ; U²+T-V )          where

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

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

-I've split this integral into 2 integrals:

   I =  §_a^b  ( y+ m ) [ ( y+m ) ( y+n ) ( y²+p.y+q ) ] ^-1/2 dy  -  §_a^b  m [ ( y+m ) ( y+n ) ( y²+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)

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 = D₀ , R03 = D_L , R04 = V_R , R05 = t₀ , R06 = R₀ ( current scale factor )

Flags:  F00  set flag F00  if you don't want to re-calculate the age of the Universe t₀

              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    D₀ = 3.401195936  E10  l-y     = R02
                                                    RDN    D_L = 4.117670678  E11  l-y     = R03
                                                    RDN    V_R = 2.482624771                   = R04         ( unit c = 1 )
                                                  LASTX   t₀  = 1.550654988  E10  years = R05

   and   R00 = k = -1  ( hyperbolic space )
           R06 =  R₀ = 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    D₀ = 2.860772072  E10  l-y     = R02
                                                    RDN    D_L = 2.283549897  E11  l-y     = R03
                                                    RDN    V_R = 2.088154797                   = R04         ( unit c = 1 )
                                                  LASTX   t₀  = 1.359670612  E10  years = R05

   and   R00 = k = +1  ( spherical space )
           R06 =  R₀ = 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