hp41programs

General Relativity and Cosmology for the HP-41

General Relativity and Cosmology for the HP-41

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

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 numbers are linked by the relation:  ( 1/H0 ) ( years ) = 977.79222143 109 / H0 ( km/s/Mpc )
-These decimals are only given to facilitate comparisons: Hubble Constant is known very inaccurately!

-The mean density of matter is of the order of  ( rho )mat ~  3  10 -28 kg/m3 , perhaps 10 or 100 times more!
( possible dark matter , non-baryonic or exotic matter ... )

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

1°) The Classical Tests

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

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:     Mercury  M = Mass of the Sun = 1  ,  T = 87.969 days  ,  e = 0.2056

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

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

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

( 30 bytes / SIZE 000 )

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

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

-And if we know  M , a , e:

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

( 30 bytes / SIZE 000 )

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

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

b) Light Deflection

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

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

( 17 bytes / SIZE 000 )

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

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

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

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

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

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

c) Gravitational Red-Shift

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

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

( 45 bytes / SIZE 000 )

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

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

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

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

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

-If both GM/(c2r) are small (  GM/(c2r) << 1 )  we can use a much shorter routine:
-Replace line 08 by 74243 E-32   if   M is expressed in kg  and  r in meter

 01  LBL "GRZ2"  02  1/X  03  CHS  04  X<>Y  05  1/X  06  +  07  *  08  2.1208 E-6  09  *  10  END

2°) Gravitational Radiation on Circular Orbits

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

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

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

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

( 32 bytes / SIZE 000 )

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

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

Example:   Let's take once again the binary pulsar PSR J0737-3039      m1 = 1.25  ,   m2 = 1.34    ( we neglect the small eccentricity )

-Currently,  r is approximately  0.00586 AU
-Find in how many time we'll have r = 0.001 AU

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

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

3°) Cosmology

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

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

ds2 = c2dt2 - R2(t)/(1-kr2/4)2 . ( dr2 + r2d(theta)2 + r2sin2(theta)d(phi)2 )

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

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

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

_________________________________________________________________________________________________________________________

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

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

_________________________________________________________________________________________________________________________

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

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

-The equations reduce to:

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

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

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

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

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

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

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

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

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

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

-Here is a summary of the different models:

********************************* MANY UNIVERSES***************************************

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

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

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

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

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

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

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

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

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

-Critical parameters:

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

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

_________________________________________________________________________________________________________________________

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

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

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

Formulae:

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

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

( 53 bytes / SIZE 000 )

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

Example:      z = 7

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

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

-Here,  k = +1

Formulae:

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

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

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

( 121 bytes / SIZE 000 )

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

Example:      q0 = 2  ;  z = 3

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

Note:

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

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

-Here,  k = -1

Formulae:

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

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

( 130 bytes / SIZE 000 )

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

Example:      q0 = 0.4  ;  z = 7

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

b) Empty Universes:   q0 + L0 = 0

b-1)  Milne Model:   q0 = 0

k = -1  ( hyperbolic geometry )

Formulae:

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

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

( 39 bytes / SIZE 000 )

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

Example:        z = 7

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

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

k = 0  ( Euclidean geometry )

Formulae:

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

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

( 39 bytes / SIZE 000 )

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

Example:        z = 7

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

b-3)  -1  <  q0  <  0

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

Formulae:       k = -1   ( hyperbolic geometry )

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

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

( 112 bytes / SIZE 000 )

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

Example:      q0 = -0.3 ;  z = 7

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

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

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

-0.3  ENTER^
7   XEQ "VAC"   >>>>  D  =  13.41113179 Gly
RDN   D0 =  32.72273316 Gly
RDN   DL = 471.6710452 Gly
RDN   VR =  2.376081547
LASTX  t0  =  15.46638674 Gy

b-4)  q0  >  0

k = -1   ( hyperbolic geometry )

Formulae:

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

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

( 101 bytes / SIZE 000 )

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

Example:      q0 =  0.3 ;  z = 7

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

b-5)  q0  <  -1

k = +1  ( spherical geometry )

Formulae:

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

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

( 107 bytes / SIZE 000 )

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

Example:      q0 =  -1.1 ;  z = 2

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

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

t0 = the time since the minimum of the scale factor R(t)
The age of the Universe is actually infinite.

Notes:

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

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

c) Eddington-Lemaitre Universes ( Critical parameters )

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

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

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

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

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

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

Data Registers: /
Flag:  F10
Subroutines: /

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

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

( 219 bytes / SIZE 000 )

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

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

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

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

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

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

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

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

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

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

d) Euclidean Universes

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

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

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

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

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

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

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

( 71 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y q0 t0 X z D

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

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

d-2)  q0 > 1/2

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

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

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

( 76 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y q0 t0 X z D

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

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

e) More General Models

e-1)  Numerical Integration ( program#1 )

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

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

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

when the program stops:

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

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

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

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

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

( 212 bytes / SIZE 017 )

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

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

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

FIX 4

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

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

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

Notes:

-If   ( 3L+2q-1)3 >= 27 L(q+L)2   and  q <= -1  there is no Big Bang  and  t0 = infinity ( 9.999999999 E99 for the HP-41 )
-Lines 86 to 103  use Gauss Legendre 3-point formula
the number n of subintervals ( in register R03 ) is doubled until 2 consecutive rounded values are equal
thus, the accuracy is controlled by the display format

-Actually, for our Universe ( q0 of the order of  -0.5 and L0 of the order of 0.6 ),
the results have at least 6 significant digits if n = 4 and lines 87 to 103 may be replaced by  4  STO 03  XEQ "GL3"  RTN

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

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

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

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

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

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

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

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

( 602 bytes / SIZE 023 )

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

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

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

-We also have

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

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

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

-We also have

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

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

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

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

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

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

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

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

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

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

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

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

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

Data Registers:   R00 thru R21

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

when the program stops:

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

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

Subroutines:

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

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

 01  LBL "COSMO3"   02  1   03  +   04  STO 16   05  FS? 00   06  GTO 04   07  CLX   08  +   09  +    10  STO 17   11  ST+ X   12  STO Z   13  +   14  1   15  -   16  STO 18   17  0   18  STO 21   19  R^   20  CHS   21  XEQ "P3"   22  STO 13   23  RDN   24  STO 14   25  X<>Y   26  STO 15   27  X<>Y   28  R^   29  1   30  XEQ 02   31  STO 19 32  90    33  TAN   34  STO 20                35  SIGN   36  CHS   37  FS? 01   38  GTO 00   39  RCL 15   40  X>Y?   41  GTO 01   42  CLX   43  RCL 14   44  X>Y?   45  GTO 01   46  X<>Y   47  LBL 00   48  RCL 13   49  XY 100  - 101  STO 10 102  RCL 15 103  RCL 14 104  RCL 13 105  RCL 16 106  XEQ 03 107  STO 12 108  RCL 15 109  RCL 14 110  RCL 13 111  1 112  XEQ 03 113  RCL 12 114  - 115  STO 02 116  RCL 10 117  STO 01 118  RCL 19 119  STO 05 120  RCL 20 121  STO 07 122  RCL 21 123  STO 08 124  2 125  RCL 17              126  / 127  SQRT 128  ST* 02 129  3 130  / 131  ST* 01 132  ST* 05 133  ST* 08 134  RCL 18 135  X#0? 136  SIGN 137  STO 00 138  RCL 02 139  STO 09 140  RCL 18 141  ABS 142  SQRT 143  X=0? 144  SIGN 145  1/X 146  STO 06 147  SF 24 148  ST* 07 149  / 150  ENTER^ 151  E^X-1 152  LASTX 153  CHS 154  E^X-1 155  - 156  2 157  / 158  RCL Y              159  R-D 160  SIN 161  RCL 00  162  X#0? 163  ENTER^ 164  X>0? 165  ENTER^ 166  R^ 167  RCL 06 168  * 169  RCL 16 170  * 171  STO 03 172  8 173  137 E8 174  LBL 05 175  ST* IND Y 176  DSE Y 177  GTO 05 178  RCL 05 179  SIGN 180  RCL 09 181  STO 04 182  RCL 03 183  RCL 02 184  RCL 01 185  CF 24 186  END

( 286 bytes / SIZE 022 )

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

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

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

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

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

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

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

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

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

f) Non-Vanishing Pressure

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

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

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

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

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

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

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

Formula:

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

( replace z by +infinity to get  t0 )

Data Registers: /
Flags: /
Subroutines: /

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

( 109 bytes / SIZE 000 )

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

( execution time = 3 seconds )

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

49 E-6  ENTER^
0.61    ENTER^
7      XEQ "TOL"  >>>>  D = 1.564459530 1010   light-years
RDN   t0 = 1.814132844 1010   years
RDN  R0 = 2.193893533 1010   light-years
RDN   k  = -1  ( hyperbolic Universe )

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

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

R10 =  z + 1                  R13 =  (Omega)rad
R11 =  (Omega)mat       R14  = (Omega)tot - 1
R12 =  (Omega)lambda  R15-R16-R17:  temp

-when the program stops:

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

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

F03  set flag F03 to fix the number n of subintervals used by "GL3"  and  store n into R09
if F03 is clear, n is doubled until 2 consecutive rounded values are equal.
F10

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

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

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

( 186 bytes / SIZE 018 )

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

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

•   for instance,  FIX 4    CF 03

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

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

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

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

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

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

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

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

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

-There are many subcases and the following program only works if the cosmological constant is positive

and if the quartic:   (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad

has 2 distinct non-positive roots and a pair of complex conjugate roots.

-These conditions are precisely satisfied by our Universe !

Formulae:

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

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

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

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

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

Data Registers:   R00 thru R32

R10 =  z + 1                  R13 =  (Omega)rad
R11 =  (Omega)mat       R14  = (Omega)tot - 1
R12 =  (Omega)lambda  R15 thru R32:  temp

-when the program stops:

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

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

F01 & F02 are used by "P4"

Subroutines:

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

-Line 124 is a three-byte  GTO 01

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

( 416 bytes / SIZE 033 )

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

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

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

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

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

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

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

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

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

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

Notes:

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

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

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

References:

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