The Museum of HP Calculators

# General Relativity and Cosmology for the HP-41

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

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

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

Overview

1°)  The Classical Tests

b)   Light Deflection
c)   Gravitational Red-shift

2°)  Gravitational Radiation on circular orbits

3°)  Cosmology

a)  Friedmann Universes:  Cosmological Constant = 0

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

b)  Empty Universes:   L0 + q0 = 0

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

c) Eddington-Lemaitre Universes ( critical parameters )

d) Euclidean Universes

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

e)  More general models

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

f)  Non-vanishing Pressure   ( new )

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

Astronomical Data:

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

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

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

1°) The Classical Tests

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

Third Kepler's law:  a3/T2 = GM/(4.PI2)

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

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

( 33 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z M / Y T M X e adv(p)

where    M = Mass ( expressed in Solar Masses )
T = Orbital Revolution ( in days )
e = eccentricity

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

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

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

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

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

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

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

( 30 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z a a Y T a X e adv(p)

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

-And if we know  M , a , e:

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

( 30 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z M / Y a / X e adv(p)

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

b) Light Deflection

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

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

( 17 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y M / X d LD(")

where  M is expressed in Solar Masses

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

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

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

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

c) Gravitational Red-Shift

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

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

( 45 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z M a Y rS a X rO adv(p)

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

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

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

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

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

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

2°) Gravitational Radiation on Circular Orbits

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

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

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

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

( 32 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T m1 m1 Z m2 m1 Y r0 m1 X r1 T

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

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

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

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

3°) Cosmology

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

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

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

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

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

_________________________________________________________________________________________________________________________

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

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

_________________________________________________________________________________________________________________________

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

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

-The equations reduce to:

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

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

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

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

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

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

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

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

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

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

-Here is a summary of the different models:

_______________________________________________MANY UNIVERSES___________________________________________________

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

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

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

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

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

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

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

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

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

-Critical parameters:

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

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

_________________________________________________________________________________________________________________________

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

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

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

Formulae:

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

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

( 53 bytes / SIZE 000 )

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

Example:      z = 7

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

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

-Here,  k = +1

Formulae:

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

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

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

( 121 bytes / SIZE 000 )

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

Example:      q0 = 2  ;  z = 3

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

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

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

-Here,  k = -1

Formulae:

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

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

( 130 bytes / SIZE 000 )

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

Example:      q0 = 0.4  ;  z = 7

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

b) Empty Universes:   q0 + L0 = 0

b-1)  Milne Model:   q0 = 0

k = -1  ( hyperbolic geometry )

Formulae:

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

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

( 39 bytes / SIZE 000 )

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

Example:        z = 7

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

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

k = 0  ( Euclidean geometry )

Formulae:

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

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

( 39 bytes / SIZE 000 )

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

Example:        z = 7

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

b-3)  -1  <  q0  <  0

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

Formulae:       k = -1   ( hyperbolic geometry )

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

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

( 112 bytes / SIZE 000 )

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

Example:      q0 = -0.3 ;  z = 7

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

b-4)  q0  >  0

k = -1   ( hyperbolic geometry )

Formulae:

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

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

( 101 bytes / SIZE 000 )

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

Example:      q0 =  0.3 ;  z = 7

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

b-5)  q0  <  -1

k = +1  ( spherical geometry )

Formulae:

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

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

( 107 bytes / SIZE 000 )

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

Example:      q0 =  -1.1 ;  z = 2

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

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

Notes:

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

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

c) Eddington-Lemaitre Universes ( Critical parameters )

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

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

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

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

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

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

Data Registers: /
Flag:  F10
Subroutines: /

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

( 219 bytes / SIZE 000 )

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

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

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

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

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

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

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

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

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

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

d) Euclidean Universes

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

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

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

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

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

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

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

( 71 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y q0 t0 X z D

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

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

d-2)  q0 > 1/2

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

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

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

( 76 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y q0 t0 X z D

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

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

e) More General Models

e-1)  Numerical Integration ( program#1 )

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

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

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

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

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

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

( 212 bytes / SIZE 017 )

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

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

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

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

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

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

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

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

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

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

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

For y = 3    3  XEQ "Y-H0T"  >>>>  2.267  and similarly:

 y 0 0.1 0.25 yI 0.75 1 1.25 1.5 1.75 2 2.5 3 10 H0.t 0 0.046 0.178 0.495 0.765 1.034 1.266 1.466 1.64 1.793 2.053 2.267 3.7

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                      ( theoretically a three-byte GTO )
132  LBL 02
133  "S2"
134  ASTO 00
135  90
136  STO 02
137  RCL 03
138  RCL 04
139  X=Y?
140  GTO 01
141  FC? 00
142  XEQ 05
143  ST+ X
144  FC? 00
145  STO 18
146  LBL 01
147  RCL 16
148  SQRT
149  1/X
150  ASIN
151  STO 02
152  FC? 00
153  XEQ 05
154  FC? 00
155  STO 17
156  RCL 10
157  RCL 16
158  *
159  SQRT
160  1/X
161  ASIN
162  GTO 04
163  LBL 03
164  RCL 11
165  1
166  +
167  FS? 01
168  X>0?
169  GTO 03
170  20
171  1/X
172  RCL 22
173  ABS
174  STO 22
175  X>Y?
176  GTO 03
177  "S3"
178  ASTO 00
179  RCL 21
180  X<>Y
181  /
182  ENTER^
183  X^2
184  1
185  +
186  SQRT
187  +
188  LN
189  CHS
190  STO 01
191  1
192  RCL 21
193  -
194  RCL 22
195  /
196  ENTER^
197  X^2
198  1
199  +
200  SQRT
201  +
202  LN
203  STO 02
204  FC? 00
205  XEQ 05
206  FC? 00
207  STO 17
208  RCL 10
209  1/X
210  RCL 21
211  -
212  RCL 22
213  /
214  ENTER^
215  X^2
216  1
217  +
218  SQRT
219  +
220  LN
221  GTO 04
222  LBL 03
223  "S"
224  ASTO 00
225  RCL 11
226  1
227  +
228  RCL 12
229  X=Y?
230  X#0?
231  GTO 00
232  90
233  TAN
234  GTO 01
235  LBL 00
236  FC? 00
237  XEQ 05
238  LBL 01
239  FC? 00
240  STO 17
241  RCL 10
242  SQRT
243  1/X
244  LBL 04
245  STO 01
246  XEQ 05
247  STO 19
248  CF 10
249  XEQ 05
250  STO 02
251  STO 20
252  RCL 14
253  ABS
254  SQRT
255  X=0?
256  SIGN
257  1/X
258  STO 03
259  STO 06
260  STO 07
261  STO 08
262  /
263  ENTER^
264  E^X-1
265  LASTX
266  CHS
267  E^X-1
268  -
269  2
270  /
271  RCL Y
272  R-D
273  SIN
274  RCL 14
275  X#0?
276  ENTER^
277  X<0?
278  ENTER^
279  GTO 07                      ( theoretically a three-byte GTO )
280  LBL 05
281   E99
282  STO 20
283  SIGN
284  STO 03
285  LBL 06
286  RCL 03
287  ST+ 03
288  XEQ "GL3"
289  VIEW X
290  RND
291  LASTX
292  X<> 20
293  RND
294  X#Y?
295  GTO 06
296  RCL 04
297  RTN
298  LBL "S"
299  X^2
300  ENTER^
301  ENTER^
302  X^2
303  RCL 12
304  *
305  RCL 14
306  +
307  *
308  RCL 13
309  +
310  SQRT                        5 bytes may be saved if you add    LBL 00   after line  336
311  1/X                            and if you replace lines 310 to 315 by  GTO 00
312  ST+ X
313  FS? 10
314  *
315  RTN
316  LBL "S1"
317  E^X
318  ENTER^
319  1/X
320  +
321  2
322  /
323  X^2
324  RCL 15
325  *
326  ENTER^
327  STO Z
328  RCL 15
329  +
330  *
331  RCL 12
332  *
333  RCL 13
334  RCL 15
335  /
336  -
337  SQRT
338  1/X
339  ST+ X
340  FS? 10
341  *
342  RTN
343  LBL "S2"
344  SIN
345  X^2
346  RCL 16
347  *
348  ENTER^
349  STO Z
350  RCL 16
351  +
352  *
353  RCL 12
354  CHS
355  *
356  RCL 13
357  RCL 16
358  /
359  +
360  SQRT
361  1/X
362  ST+ X
363  D-R
364  FS? 10
365  *
366  RTN
367  LBL "S3"
368  E^X
369  ENTER^
370  1/X
371  -
372  2
373  /
374  RCL 22
375  *
376  RCL 21
377  +
378  STO Y
379  LASTX
380  ST+ X
381  +
382  RCL 12
383  *
384  FS? 10
385  GTO 00
386  *
387  SQRT
388  1/X
389  RTN
390  LBL 00
391  /
392  SQRT
393  RTN
394  LBL 07
395  R^
396  RCL 10
397  *
398  ST* 03
399  RCL 14
400  X#0?
401  SIGN
402  CHS
403  STO 00
404  RCL 17
405  STO 05
406  RCL 19
407  STO 01
408  RCL 15
409  ST* 07
410  SF 24
411  RCL 16
412  ST* 08
413  RCL 18
414  STO 09
415  9
416  137 E8
417  LBL 08
418  ST* IND Y
419  DSE Y
420  GTO 08
421  RCL 05
422  SIGN
423  RCL 20
424  STO 04
425  RCL 03
426  RCL 02
427  RCL 01
428  CLA
429  CLD
430  CF 24
431  END

( 602 bytes / SIZE 023 )

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Data Registers:   R00 thru R21

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

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

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

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

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

( 286 bytes / SIZE 022 )

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

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

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

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

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

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

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

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

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

f) Non-Vanishing Pressure

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

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

where          a = Stefan's Constant  and  T = temperature of the Cosmic Microwave Background

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

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

-We now have  (Omega)mat = 2(q0+L00)

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

Formula:

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

( replace z by +infinity to get  t0 )

Data Registers: /
Flags: /
Subroutines: /

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

( 111 bytes / SIZE 000 )

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

( execution time = 3 seconds )

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

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

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

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

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

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

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

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

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

( 200 bytes / SIZE 019 )

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

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

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

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

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

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

-These conditions are precisely satisfied by our Universe!

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

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

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

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

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

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

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

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

Data Registers:   R00 thru R34

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

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

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

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

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

( 424 bytes / SIZE 035 )

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

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

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

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

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

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