Template

# Age of the Universe & Redshift -> Distance for the HP-41

Overview

1°) Empty Universes

a)  0 < lambda < 1
b)  General Case

2°) Tolman Universes

a)  General Case
b)  Cyclic Universes

3°) Our Universe ?

a)  Via Carlson Elliptic Integrals ( 2 programs )
b)  Via Numerical Integration
( 3 programs )

4°)  More General Cyclic Universes  ( 3 programs )

5°)  More General Universes

a)  Cosmological Constant = 0  ( 2 programs )
b)  Cosmological Constant # 0
( 4 programs )

6°)  Redshifht -> Distance

a)  Gauss-Legendre 2-point formula
c)  Cosmological Constant = 0 , Positive Matter, Negative Pressure
d)  No Matter, Negative Pressure
e)  No Pressure, Negative Matter

Latest Update:   paragraph 6°)

-The age of the universe may be computed by the following integral:

t = (1/H)       § 0 1  y. [  (Omega) lambda.y 4   + ( 1-(Omega)tot ).y2 +  (Omega) mat .y + (Omega)rad  ] -1/2  dy

where  (Omega)tot = (Omega)mat + (Omega) lambda + (Omega)rad     ,   H = Hubble constant

-When the minimum radius of the Universe is > 0, the age of the Universe is infinite, but these programs calculate the time since the last minimum.

-In several programs below  1/H = 13.77172143  E9  years  which corresponds to H = 71 km/s/Mpc
-More generally,  1/H ( in GigaYears ) = 977.7922214 / H  ( if H is expressed in km/s/Mpc )

-In paragraph 6, we calculate:

D   =  light-travel time distance ( in Giga-light-years ) = light-travel time ( in Giga-years )
D0 = comoving distance ( in Giga-light-years )
DL = luminosity-distance ( in Giga-light-years )

Formulae:

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 = 1/yem  ,  z = redshift  ,   y = R/R0

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

F(x) = Sinh(x)  if  k = -1     hyperbolic space
DL = R0 ( z + 1 ) F(D0/R0)    where    F(x) =   x         if  k = 0      euclidean space                      R0 = current scale factor ( or radius )
F(x) = Sin(x)    if  k = +1    spherical space

1°) Empty Universes

a)  0 < Lambda < 1

-Here,  (Omega)mat = (Omega)rad = 0

Data Registers: /

Flags: /
Subroutines: /

 01 LBL "AEU" 02 SQRT 03 STO Y 04 1 05 LASTX 06 - 07 SQRT 08 1/X 09 STO T 10 * 11 ENTER 12 X^2 13 SIGN 14 ST+ L 15 X<> L 16 SQRT 17 + 18 LN 19 X<>Y 20 / 21 13.77172143 E9 22 ST* Z 23 * 24 END

( 49 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y / R X (Omega)lambda t

Example:
(Omega)lambda = 0.41

0.41   XEQ "AEU"   >>>>   16318046180  years
X<>Y   17929254160  light-years

Note:

-The constant k is always equal to  -1  in these cases ( hyperbolic universes )

b)  General Case

Data Registers: /

Flags: /
Subroutines: /

 01 LBL "AEU" 02 DEG 03 X#0? 04 GTO 00 05 CLST 06 SIGN 07 CHS 08 13.77172143 E9 09 ENTER 10 RTN 11 LBL 00 12 ENTER 13 CF 06 14 X<0? 15 SF 06 16 ABS 17 SQRT 18 1/X 19 LASTX 20 1 21 R^ 22 X#Y? 23 GTO 00 24 CLST 25 13.77172143 E9 26 90 27 TAN 28 R^ 29 + 30 RTN 31 LBL 00 32 - 33 CF 07 34 X<0? 35 SF 07 36 ABS 37 SQRT 38 1/X 39 STO Z 40 * 41 FS? 06 42 GTO 01              43 ENTER 44 X^2 45 SIGN 46 FS? 07 47 CHS 48 ST+ L 49 X<> L 50 SQRT 51 + 52 LN 53 GTO 02 54 LBL 01 55 ASIN 56 D-R 57 LBL 02 58 R^ 59 * 60 13.77172143 E9 61 ST* Z 62 ST* T 63 * 64 R^ 65 SIGN 66 CLX 67 PI 68 FC? 06 69 CLX 70 ST* L 71 CLX 72 FC?C 06 73 FS? 07 74 X<0? 75 STO T               76 ISG X 77 CLX 78 FC?C 07 79 CHS 80 X<> Z 81 X<>Y 82 END

( 151 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T / Rmin or  R max Z / k ( -1 , 0 , +1 ) Y / R X (Omega)lambda t L / 0 or P

Where  P  is the period

Example1:

0.4   CHS   XEQ "AEU"    >>>>>     t = 1.2279852 E10  years
RDN      R = 1.1639228 E10  light-years
RDN       k = -1
RDN  Rmax = 2.1775003 E10  light-years
LASTX    P = 6.8408191 E10  years

Example2:

0.4   XEQ "AEU"    >>>>>     t = 1.6233224 E10  years
RDN      R = 1.7779215 E10  light-years
RDN       k = -1
RDN  Rmin =  0
LASTX    P =  0   ( no period )

Example3:

1.4   XEQ "AEU"    >>>>>     t = 1.4420357 E10  years
RDN      R = 2.1775003 E10  light-years
RDN       k = +1
RDN  Rmin = 1.1639228 E10  light-years
LASTX    P =  0   ( no period )

Example4:

0   XEQ "AEU"    >>>>>     t = 1.3771721 E10  years
RDN     R = 1.3771721 E10  light-years
RDN      k = -1
RDN  Rmin = 0
LASTX    P =  0   ( no period )

Example5:

1   XEQ "AEU"    >>>>>     t = 9.9999999 E99  years  ( in fact: infinite )
RDN     R = 1.3771721 E10  light-years
RDN      k = 0
RDN  Rmin = 0
LASTX   P =  0   ( no period )

2°) Tolman Universes

a)  General Case

-Here,  (Omega)mat = 0
-Such models have been studied by Tolman.

Data Registers:    R00 = k ( -1 , 0 , +1 )               R06 to R10:  temp

R01 = t            R03 = Rmin       R05 = Rmax
R02 = P           R04 = R
Flags:  F24
Subroutines:

 01 LBL "AUR" 02 DEG 03 STO 06 04 X<>Y 05 STO 07 06 + 07 90 08 TAN 09 STO 02 10 STO 05 11 STO 09              12 CLX 13 STO 03 14 SIGN 15 - 16 STO 00 17 CHS 18 RCL 06 19 X=0? 20 SIGN 21 ST/ 07 22 / 23 STO 08 24 RCL 06 25 X=0? 26 GTO 01 27 ABS 28 SQRT 29 ST+ X 30 STO 10 31 CLX 32 2 33 / 34 STO 08 35 X^2 36 RCL 07 37 - 38 X<0? 39 GTO 02 40 SQRT 41 RCL 08 42 SIGN 43 * 44 RCL 08              45 + 46 CHS 47 RCL 07 48 RCL Y 49 X=0? 50 SIGN 51 / 52 X<0? 53 CLX 54 SQRT 55 X<>Y 56 X<0? 57 CLX 58 SQRT 59 XY 61 GTO 03 62 LBL 01 63 RCL 07 64 CHS 65 RCL 08 66 X=0? 67 GTO 02 68 / 69 X<0? 70 CLX 71 SQRT 72 ENTER 73 GTO 03 74 LBL 02 75 CLST 76 LBL 03 77 1 78 X<>Y 79 X>Y? 80 STO 05 81 X<> Z 82 X>Y? 83 STO 05              84 X>0? 85 X>Y? 86 FS? 30 87 STO 03 88 X<> Z 89 X>0? 90 X>Y? 91 FS? 30 92 STO 03 93 RCL 06 94 X=0? 95 GTO 01 96 X<0? 97 GTO 02 98 RCL 07 99 RCL 08 100 ST+ X 101 + 102 1 103 + 104 SQRT 105 1 106 + 107 RCL 08 108 ST+ Y 109 RCL 07 110 X<0? 111 CLX 112 SQRT 113 + 114 X=0? 115 GTO 00 116 STO 04              117 / 118 ABS 119 LN 120 RCL 10 121 / 122 STO 01 123 RCL 09 124 RCL 05 125 X=Y? 126 GTO 03 127 X^2 128 RCL 08 129 + 130 RCL 04 131 / 132 ABS 133 LN 134 ST+ X 135 RCL 10 136 / 137 STO 02 138 GTO 03 139 LBL 01 140 1 141 RCL 07 142 X<0? 143 CLX 144 SQRT 145 STO 04 146 - 147 RCL 08 148 X=0? 149 GTO 04 150 / 151 STO 01              152 RCL 09 153 RCL 05 154 X=Y? 155 GTO 03 156 X^2 157 RCL 08 158 * 159 RCL 07 160 + 161 SQRT 162 RCL 04 163 - 164 ST+ X 165 RCL 08 166 / 167 STO 02 168 GTO 03 169 LBL 02 170 1 171 RCL 08 172 ST+ Y 173 X^2 174 RCL 07 175 - 176 SQRT 177 / 178 ASIN 179 90 180 + 181 D-R 182 RCL 10              183 / 184 STO 01 185 PI 186 ST+ X 187 LASTX 188 / 189 STO 02 190 GTO 03 191 LBL 04 192 .5 193 STO 01 194 GTO 03 195 LBL 00 196 RCL 09 197 STO 01 198 LBL 03 199 RCL 00 200 X#0? 201 SIGN 202 ENTER 203 X<> 00 204 ABS 205 SQRT 206 X=0? 207 SIGN 208 1/X 209 SF 24 210 13.77172143 E9 211 ST* 01 212 ST* 02 213 * 214 STO 04              215 ST* 03 216 ST* 05 217 RCL 01 218 RCL 02 219 RDN 220 CF 24 221 END

( 279 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS T / P Z / k Y (Omega)rad R X (Omega)lambda t

Example1:         (Omega)lambda   = 0.732 ,   (Omega)rad  = 0. 000049

0.000049   ENTER^
0.732         XEQ "AUR"   >>>>    t  = 20197856760 years         = R10
RDN    R = 26604833060 light-years = R04
RDN    k = -1                                     = R00
RDN    P = 9.999999999 E99            = R02 = period  ( infinity = no period )

-We also have   R03 = Rmin =  0   &   R05 = Rmax = 9.999999999 E99

Example2:         (Omega)lambda   = -0.1 ,   (Omega)rad  = -0.001

0.001  CHS   ENTER^
0.1      CHS    XEQ "AUR"   >>>>    t  = 13327008450 years        = R10
RDN    R = 13124856690 light-years = R04
RDN    k = -1                                     = R00
RDN    P = 136816382100 years       = R02 = period

-We also have   R03 = Rmin =  395565882.1 light-years   &   R05 = Rmax = 43548210540  light-years

Notes:

-The constant k = -1 , 0 , +1  corresponds to hyperbolic, euclidean, spherical universes respectively.
-This program also works for pulsating universes.

b)  Cyclic Universes

-With this program, you choose the minimum, current and maximum scale factor:  Rmin R0 Rmax
and "AUR" returns the age of the Universe T, the period P, the current Hubble "constant" H0 and the deceleration parameter q

Data Registers:       R00 = L (Gy)-2

R01 = Rmin      R04 = T       R07 = q
R02 = R0         R05 = P
R03 = Rmax     R06 = H0 ( km/s/Mpc )
Flags: /
Subroutines: /

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

( 88 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS T / q Z Rmin H0 Y R0 P X Rmax T

Where  T = age of the Universe ,  P  = the period ,  H0 = the current Hubble "constant" ( in km/s/Mpc ) and  q = the deceleration parameter

Example:
Rmin = 1    R0 = 14   Rmax = 314   ( Gygalightyears )

1    ENTER^
14   ENTER^
314  XEQ "AUR"  >>>>   T  = 13.96898880 = R04               ---Execution time = 3s---
RDN   P  = 986.4650960 = R05
RDN  H0 = 69.59427440 = R06
RDN   q = -0.003136335 = R07

-We also have the cosmological constant in  R00 = L = -0.000030427 (Gy)-2

Note:

k is always equal to -1  ( hyperbolic space )

3°) Our Universe ?

a) Via Carlson Integrals

-The following program only works if the cosmological constant is positive

and if the quartic:   (Omega)lambda.y 4 + ( 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 !

Data Registers:     R00 thru R21: temp
Flags:  F00-F01-F02-F10
Subroutines:  "RFZ"  "RJZ"  "P4"   (  cf "Elliptic Integrals for the HP41" & "Polynomials for the HP41" )

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

( 230 bytes / SIZE 022 )

 STACK INPUTS OUTPUTS Z (Omega) rad k Y (Omega)lamda R X (Omega) mat t

Example1:      (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega)rad  = 0.000049

0.000049   ENTER^
0.521         ENTER^
0.044         XEQ "AUM"    >>>>   t = 15602214530 years
RDN   R= 20881806280 light-years
RDN   k = -1

Example2:      (Omega)mat = 0 .271 ,  (Omega)lambda  = 0.732 ,  (Omega) rad  = 0.000049

0.000049   ENTER^
0.732         ENTER^
0.271         XEQ "AUM"    >>>>   t =  13667886770 years
RDN   R= 249407504600 light-years
RDN   k = +1

Note:

-We can also use M-Code routines  RF & RJ
-The following program employs a shorter formula given in reference [1]

Data Registers:     R00 thru R15: temp
Flags:  F00-F01-F02-F10
Subroutines:  "RFZ"  "RJZ" ( M-Code ) and  "P4"   (  cf "Carlson Elliptic Integrals" & "Polynomials for the HP41" )

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

( 203 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T H0 / Z (Omega) mat k Y (Omega)lamda R X (Omega)rad t

Example1:    H0 = 71 km/s/Mpc     (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega)rad  = 0.000049

71           ENTER^
0.044         ENTER^
0.521         ENTER^
0.000049      XEQ "AUM"    >>>>   t = 15.60221456 Giga-years                                         ---Execution time = 34s---
RDN   R= 20.88180628 Giga-light-years
RDN   k = -1

Example2:     H0 = 71 km/s/Mpc      (Omega)mat = 0 .271 ,  (Omega)lambda  = 0.732 ,  (Omega) rad  = 0.000049

71        ENTER^
0.271      ENTER^
0.732      ENTER^
0.000049  XEQ "AUM"    >>>>   t = 13.66788679  Giga-years                                         ---Execution time = 34s---
RDN   R= 249.4075046 Giga-light-years
RDN   k = +1

b) Via Numerical Integration  ( 3 programs )

-This routine employs Gauss-Legendre 3-point formula to evaluate the integral.
-The formula is applied to N subintervals

Data Registers:     R00 & R04 thru R10: temp                                  ( Registers R01 thru R03 are to be initialized before executing "AUM" )

•  R01 = (Omega)mat
•  R02 = (Omega)lambda
Flags: /
Subroutines: /

 01 LBL "AUM" 02 STO 05 03 1/X 04 STO 06       05 2 06 / 07 STO 07 08 1.6 09 STO 10 10 FRC 11 SQRT 12 * 13 STO 08 14 RCL 01 15 RCL 02 16 + 17 RCL 03       18 + 19 1 20 - 21 STO 00 22 CLX 23 STO 04 24 GTO 01 25 LBL 00 26 STO 09 27 X^2 28 ENTER 29 ENTER 30 X^2 31 RCL 02 32 * 33 RCL 00       34 - 35 * 36 RCL 01 37 + 38 * 39 RCL 03 40 + 41 SQRT 42 / 43 RCL 09 44 * 45 RTN 46 LBL 01 47 RCL 07 48 RCL 08       49 - 50 XEQ 00 51 ST+ 04 52 RCL 07 53 XEQ 00 54 RCL 10 55 * 56 ST+ 04 57 RCL 07 58 RCL 08 59 + 60 XEQ 00 61 ST+ 04 62 RCL 06       63 ST+ 07 64 DSE 05 65 GTO 01 66 RCL 00 67 X#0? 68 SIGN 69 RCL 00 70 ABS 71 SQRT 72 X=0? 73 SIGN 74 1/X 75 1.8 76 * 77 RCL 04       78 R^ 79 * 80 7650956350 81 ST* Z 82 * 83 END

( 118 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS Z / k Y / R X N t

Where N = number of subintervals

Example1:      (Omega)mat = 0 .044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049

0.000049   STO 03
0.521         STO 02
0.044         STO 01      and if you choose  N = 50

50    XEQ "AUM"       >>>>   t = 15602214540 years
RDN   R= 20881806270 light-years
RDN   k = -1

Example2:      (Omega)mat = 0 .271 ,  (Omega)lambda  = 0.732 ,  (Omega) rad  = 0.000049

0.000049   STO 03
0.732         STO 02
0.271         STO 01      and if you choose  N = 50

50         XEQ "AUM"      >>>>   t =  13667886820 years
RDN   R= 249407504600 light-years
RDN   k = +1

Notes:

-N = 50 produces a very good precision in most cases.
-However, this program is slower with N = 50 than the version listed in paragraph 3°) a)
-On the other hand, Hubble constant is not known with an accuracy of 10 digits.
-So a smaller N-value is enough to get a 4-digit precision:   N = 10  seems quite enough !

-We can also use Gauss-Legendre 2-point formula.
-The following program makes the change of variable:     y (Omega)mat + (Omega) rad  =  u2

Data Registers:     R00 = H                               R04 thru R10: temp

R01 = (Omega)mat
R02 = (Omega)lambda
Flags: /
Subroutines: /

 01 LBL "AUM"  02 STO 03  03 X<>Y  04 STO 02  05 +  06 X<>Y  07 STO 01            08 +  09 X<>Y  10 STO 00  11 SIGN  12 -  13 STO 04  14 RCL 01  15 RCL 03 16 +  17 SQRT  18 RCL 03  19 SQRT  20 STO 07  21 -  22 50  23 STO 09            24 ST+ 09    25 /  26 STO Y  27 3  28 SQRT  29 /  30 STO 08 31 -  32 STO 05  33 2  34 /  35 ST- 07  36 CLX  37 STO 06            38 LBL 01  39 RCL 07  40 RCL 08  41 X<> 05  42 STO 08  43 +  44 STO 07  45 X^2 46 STO 10  47 RCL 03  48 -  49 RCL 01  50 /  51 ENTER  52 X^2  53 RCL X  54 RCL 02            55 *  56 RCL 04  57 -  58 *  59 RCL 10  60 + 61 SQRT  62 /  63 RCL 07  64 *  65 ST+ 06  66 DSE 09  67 GTO 01  68 RCL 06  69 RCL 08  70 *  71 RCL 01            72 /  73 3  74 SQRT  75 *  76 977.7922214 77 RCL 00  78 /  79 STO 05  80 *  81 RCL 04  82 X#0?  83 SIGN  84 RCL 05  85 RCL 04            86 ABS  87 SQRT  88 X=0?  89 SIGN  90 /  91 R^  92 END

( 119 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS T H ( km/s/Mpc ) / Z (Omega)mat > 0 k Y (Omega)lambda R X (Omega)rad T

Example1:
H = 71 km/s/Mpc  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049

71       ENTER^
0.044     ENTER^
0.521     ENTER^
0.000049   XEQ "AUM"   >>>>   T = 15.60221456   ( Gy )                                         ---Execution time = 107s---
RDN    R = 20.88180628  ( Gly )
RDN    k = - 1  ( hyperbolic space )

Example2:    H = 71 km/s/Mpc   (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,  (Omega) rad  = 0.000049

71       ENTER^
0.271     ENTER^
0.732     ENTER^
0.000049   XEQ "AUM"   >>>>   T = 13.66788680   ( Gy )                                         ---Execution time = 107s---
RDN    R = 249.4075046  ( Gly )
RDN    k = + 1  ( spherical space )

Example3:    H = 71 km/s/Mpc   (Omega)mat = 0.271 ,  (Omega)lambda  = 0.729 ,  (Omega) rad  = 0

71       ENTER^
0.271     ENTER^
0.729     ENTER^
0        XEQ "AUM"   >>>>   T = 13.65710890   ( Gy )                                         ---Execution time = 107s---
RDN    R = 13.77172143  ( Gly )
RDN    k =  0  ( Euclidean space )

Notes:

50 subintervals seem enough for a very good precision
Change line 22 if you prefer another value.

-We may simplify the program if  (Omega) rad  = 0
-The following routine employs again Gauss-Legendre 2-point formula ( and the change of variable  y2 = u )

Data Registers:   R00 to R07: temp
Flags: /
Subroutines: /

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

( 107 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS Z H ( km/s/Mpc ) k Y q R X (Omega)mat T

Example1:
H = 71 km/s/Mpc   q = -0.5    (Omega)mat = 0.044

71       ENTER^
-0.5       ENTER^
0.044     XEQ "AUM"   >>>>    T  = 15.62753949   ( Gy )                                         ---Execution time = 82s---
RDN    R = 20.90467233  ( Gly )
RDN    k = - 1  ( hyperbolic space )

Example2:    H = 71 km/s/Mpc   q = -0.6   (Omega)mat = 0.271

71       ENTER^
-0.6      ENTER^
0.271    XEQ "AUM"   >>>>   T = 13.68753141   ( Gy )                                         ---Execution time = 82s---
RDN    R = 170.8171812  ( Gly )
RDN    k = + 1  ( spherical space )

Example3:    H = 71 km/s/Mpc    q = -0.595    (Omega)mat = 0.27

71      ENTER^
-0.595    ENTER^
0.27     XEQ "AUM"   >>>>   T = 13.67100708   ( Gy )                                         ---Execution time = 82s---
RDN    R = 13.77172143  ( Gly )
RDN    k =  0  ( Euclidean space )

Example4:    H = 71 km/s/Mpc    q = -0.5    (Omega)mat = 0

71     ENTER^
-0.5    ENTER^
0      XEQ "AUM"   >>>>   T = 17.16576879   ( Gy )                                         ---Execution time = 82s---
RDN    R = 19.47615522  ( Gly )
RDN    k = - 1  ( hyperbolic space )

Notes:

-N = 50 subintervals seem enough for a very good accuracy ( line 20 )

-Here we have  q = (1/2) (Omega)mat -  (Omega)lambda

-Though this program works well if  (Omega)mat = 0 , it's much slower than the routine listed in paragraph 1°)

4°)  More General Cyclic Universes

-The age of the Universe  T = (c/H)   § ymin 1  y. [  (Omega) lambda.y 4   + ( 1-(Omega)tot ).y2 +  (Omega) mat .y + (Omega)rad  ] -1/2  dy
-The Period  P = 2 (c/H) § ymin ymax  y. [  (Omega) lambda.y 4   + ( 1-(Omega)tot ).y2 +  (Omega) mat .y + (Omega)rad  ] -1/2  dy

-The following program assumes that the quartic equation    (Omega) lambda.y 4   + ( 1-(Omega)tot ).y2 +  (Omega) mat .y + (Omega)rad  = 0  has (at least ) 2 positive real roots.
-They are computed by Newton's method.

-The integrals are then calculated by Gauss-Legendre 3-point formula, applied to N subintervals between ymin & 1  and 3.N subintervals between 1 & ymax ( lines 266-267 )

Data Registers:     R00 = H ( km/s/Mpc )                                                         ( Registers R00 thru R03 are to be initialized before executing "ACU" )

•  R01 = (Omega)mat                     R04 = Rmin                        R06 thru R18: temp
•  R02 = (Omega)lambda < 0          R05 = Rmax
Flags: /
Subroutines: /

 01 LBL "ACU"  02 STO 16  03 1  04 RCL 01            05 -  06 RCL 02  07 -  08 RCL 03  09 -  10 STO 06  11 STO 17  12 2 E-8  13 STO 11  14 CLX  15 STO 09  16 STO 10  17 LBL 01  18 RCL 10  19 ENTER  20 ENTER  21 X^2  22 RCL 02  23 *  24 RCL 06  25 +  26 *  27 RCL 01  28 +  29 *  30 RCL 03            31 +  32 X<>Y  33 X^2  34 RCL 02  35 *  36 ST+ X  37 RCL 06  38 +  39 R^  40 *  41 ST+ X  42 RCL 01 43 +  44 /  45 -  46 STO 10            47 ST- Y  48 X=0?  49 SIGN  50 /  51 ABS  52 RCL 11  53 XY  83 3  84 * 85 RCL 04  86 ST+ X  87 +  88 RCL 02            89 *  90 R^  91 *  92 RCL 07  93 +  94 /  95 -  96 STO 10  97 ST- Y  98 X=0?  99 SIGN 100 / 101 ABS 102 RCL 11 103 XY? 135 X<>Y 136 RCL 04           137 X>Y? 138 X<>Y 139 RCL 05 140 X>Y? 141 X<>Y 142 STO 07 143 RDN 144 XY 146 RDN 147 XY 149 R^ 150 XY 152 STO 05 153 RDN 154 STO 04 155 X<>Y 156 STO 06 157 1 158 RCL 04 159 X 04 163 X<> 05 164 STO 06 165 LBL 03 166 RCL 06 167 RCL 07 168 + 169 CHS 170 X<> 06 171 ST* 07 172 GTO 05 173 LBL 04 174 RCL 04 175 RCL 05 176 X>Y? 177 X<>Y 178 STO 04           179 X<>Y 180 STO 05 181 + 182 STO 06 183 LBL 05 184 1 185 RCL 05 186 RCL 04 187 X<=0? 188 LN 189 ST- Z 190 - 191 STO 08 192 / 193 SQRT 194 RAD 195 ASIN 196 STO 11 197 1.6 198 STO 15 199 GTO 14 200 LBL 00 201 SIN 202 X^2 203 RCL 08 204 * 205 RCL 04 206 + 207 ENTER 208 STO Z 209 RCL 06 210 + 211 * 212 RCL 07 213 + 214 SQRT 215 / 216 RTN 217 LBL 10 218 RCL 11 219 RCL 09           220 - 221 X<>Y 222 STO 13 223 / 224 STO 14 225 2 226 / 227 ST+ 09 228 .6 229 SQRT 230 * 231 STO 10 232 CLX 233 STO 12 234 LBL 12 235 RCL 09 236 RCL 10 237 - 238 XEQ 00 239 ST+ 12 240 RCL 09 241 XEQ 00 242 RCL 15 243 * 244 ST+ 12 245 RCL 09 246 RCL 10 247 + 248 XEQ 00 249 ST+ 12 250 RCL 14 251 ST+ 09 252 DSE 13 253 GTO 12 254 RCL 12 255 * 256 RTN 257 LBL 14 258 RCL 16 259 XEQ 10 260 STO 18           261 SIGN 262 ASIN 263 X<> 11 264 STO 09 265 RCL 16 266 3 267 * 268 XEQ 10 269 RCL 18 270 ST+ Y 271 X<>Y 272 ST+ X 273 1.8 274 RCL 02 275 CHS 276 SQRT 277 * 278 ST/ Z 279 / 280 977.7922214 281 RCL 00 282 / 283 ST* Y 284 ST* Z 285 RCL 17 286 SQRT 287 / 288 ST* 04 289 ST* 05 290 X<> Z 291 DEG 292 END

( 370 bytes / SIZE 019 )

 STACK INPUTS OUTPUTS Z / R Y / P X N T

Where  N = number of subintervals   ,  T = Age of the Universe in Gigayears  ,
P = Period of the Universe in Gigayears   ,   R = radius of the Universe in Gigalightyears

Example:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001

71      STO 00

-0.1     STO 01
-0.003   STO 02
-0.001   STO 03      and if you choose  N = 8   subintervals

8       XEQ "ACU"    >>>>   T = 14.61256009   Gy                                           ---Execution time = 2mn10s---
RDN    P = 796.1455569   Gy
RDN    R =  13.10701187  GLy

-We also have:   R04 = minimum radius of the Universe  = Rmin = 1.306380222  GLy
R05 = maximum radius of the Universe = Rmax = 250.8400402 GLy

Notes:

-The numerical integration is performed with the change of argument:  y = A + ( B - A ) Sin2 u   where  A = Rmin/R  &  B = Rmax/R
-We must choose:   1-(Omega)tot > 0
-The constant k = -1  ( hyperbolic universes )

-We can also choose the minimum, current & maximum scale factors  Rmin R0 Rmax
-The following program employs Carlson elliptic integral  RJ

Data Registers:     R00 =  L = cosmological constant ( Gy-2 )

R01 = Rmin            R04 = H0                        R06 thru R15: temp
R02 = R0               R05 =  q
R03 = Rmax
Flag:  F07
Subroutines: /

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

( 288 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T / q Z Rmin H0 Y R0 P X Rmax T

Where  T = age of the Universe ,  P  = the period ,  H0 = the current Hubble "constant" and  q = the deceleration parameter

Example:
Rmin = 1    R0 = 14   Rmax = 314   ( Gigalightyears )

1    ENTER^
14   ENTER^
314  XEQ "AUC"   >>>>   T = 15.49057984   gigayears                                 ---Execution time = 58s---
RDN   P = 993.8651299   gigayears
RDN  H0 = 67.22990237  km/s/Mpc
RDN  q  =  -0.036404801

and  R00 = L = -0.000030330 Gy-2

-We can use M-Code routines  RF  &  RJ

-I've also remarked that the change of variable:    x = (a.y+r) / (y+1)   yields

§ar x dx / sqrt [ x (b-x) (x+c) (x+d) ]  =  §0+infinity  sqrt(r-a)  sqrt (a.y+r)  /  (y+1) / sqrt [ ((b-a).y + (b-r)) ((a+c).y + (r+c)) ]  dy

=  ...  =  2 SQRT [ (r-a) / a / (b-a) / (a+c) ]  {  a RF [ (b-r)/(b-a) , (r+c)/(a+c) , r/a ]  + [(r-a)/3]  RJ [ (b-r)/(b-a) , (r+c)/(a+c) , r/a ; 1 ] }

and with this formula, the program becomes slightly shorter.

Data Registers:     R00 =  L = cosmological constant ( Gy-2 )

R01 = Rmin            R04 = H0           R10 = P             R06 thru R11: temp
R02 = R0               R05 =  q
R03 = Rmax
Flag:  F07
Subroutines:  M-Code routines:  RF & RJ  ( cf "Carlson Elliptic Integrals for the HP41" )

 01 LBL "AUM"  02 STO 03            03 STO 04  04 STO 05  05 RDN  06 STO 02  07 X<>Y  08 STO 01  09 ST- 04  10 ST+ 05  11 RCL 05  12 +  13 STO 09  14 RCL 03  15 SF 07  16 XEQ 01  17 ST+ X  18 STO 10  19 RCL 02 20 LBL 01  21 STO 00            22 RCL 01  23 /  24 STO 11  25 RCL 03  26 RCL 00  27 -  28 RCL 04  29 /  30 STO 06  31 RCL 00  32 RCL 05  33 +  34 RCL 09  35 /  36 STO 07  37 RF  38 RCL 01 39 *  40 STO 08            41 SIGN  42 RCL 06  43 RCL 07  44 RCL 11  45 RJ  46 RCL 00  47 RCL 01  48 -  49 STO 06  50 *  51 3  52 /  53 RCL 08  54 +  55 RCL 06  56 RCL 01  57 / 58 RCL 04  59 /  60 RCL 09            61 /  62 SQRT  63 ST+ X  64 *  65 FS?C 07  66 RTN  67 RCL 03  68 RCL 05  69 *  70 RCL 01  71 X^2  72 +  73 STO 00  74 SQRT  75 ST* 10  76 * 77 3  78 RCL 00            79 /  80 CHS  81 STO 00  82 RCL 02  83 ST+ 05  84 RCL 01  85 -  86 *  87 RCL 02  88 RCL 03  89 -  90 *  91 RCL 05  92 *  93 RCL 02  94 3  95 ST/ Z  96 Y^X 97 /  98 STO 04            99 RCL 02 100 X^2 101 1/X 102 - 103 RCL 00 104 - 105 RCL 04 106 ST+ X 107 / 108 STO 05 109 RCL 04 110 SQRT 111 977.7922214 112 * 113 STO 04 114 RCL 10 115 R^ 116 END

( 148 bytes / SIZE 012 )

 STACK INPUTS OUTPUTS T / q Z Rmin H0 Y R0 P X Rmax T

Where  T = age of the Universe ,  P  = the period ,  H0 = the current Hubble "constant" and  q = the deceleration parameter

Example:
Rmin = 1    R0 = 14   Rmax = 314   ( Gigalightyears )

1    ENTER^
14   ENTER^
314  XEQ "AUC"   >>>>   T = 15.49057983   gigayears                                 ---Execution time = 42s---
RDN   P = 993.8651299   gigayears
RDN  H0 = 67.22990236  km/s/Mpc
RDN  q  =  -0.036404801

and  R00 = L = -0.000030330 Gy-2

Note:

-Cf also "Cyclic Universes for the HP41"

4°) More General Universes

a)  Cosmological Constant = 0

-This program uses Einstein's equation with  Lambda = 0  &  k = +1 ( spherical space )
-You choose Rmin R0 Rmax and "AUL" returns the age of the Universe T, the period P, the current Hubble "constant" H0 and the deceleration parameter q

Data Registers:       R00 = H0 ( km/s/Mpc )

R01 = Rmin      R04 = T
R02 = R0         R05 = P
R03 = Rmax     R06 = q
Flags: /
Subroutines: /

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

( 97 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS T / q Z Rmin H0 Y R0 P X Rmax T

Where  T = age of the Universe ,  P  = the period ,  H0 = the current Hubble "constant" and  q = the deceleration parameter

Example:
Rmin = 2    R0 = 41   Rmax = 257   ( Gygalightyears )

2    ENTER^
41   ENTER^
257  XEQ "AUL"   >>>>   T = 12.28422138   gigayears                                 ---Execution time = 4s---
RDN   P = 813.6724974   gigayears
RDN  H0 = 53.38731059  km/s/Mpc
RDN  q  =  0.569266382

Note:

-The following program is more general:
-You choose the parameters Hubble "constant" ,  (Omega)mat &  (Omega) rad

-Since  (Omega) lambda = 0  ,  we can calculate the integrals with elementary functions ( cf reference [2] )

Data Registers:    R00 = H0

R01 =
(Omega)mat      R03 = T          R05 = Rmin      R07 = Rmax      R09-R10-R11:  temp
R02 =
(Omega)rad       R04 = P          R06 = R0         R08 = k
Flags:  F24
Subroutines:  /

-Lines 27-78-102  are three-byte GTOs

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

( 313 bytes / SIZE 012 )

 STACK INPUTS OUTPUTS T / k Z H0 R0 Y (Omega)mat P X (Omega)rad T

Where    T = age of the Universe ( or the time since the last minimum scale factor )   P = period  ( or 0 if no period )   R0 = current scale factor

Example1:
H0 = 68 km/s/Mpc     (Omega)mat = 1.4       (Omega)rad = -0.2

68  ENTER^
1.4  ENTER^
-0.2  XEQ "AUL"   >>>>    T  = 10.22400044  Gy                                       ---Execution time = 5s---
RDN     P  = 707.0833001  Gy
RDN    R0 =  32.15308639  Gly
RDN    k  =  +1  ( spherical space )

and    R05 = 4.691072085 Gly  = minimum radius of this Universe
R07 = 220.3805326 Gly = maximum radius of this Universe.

Example2:    H0 = 68 km/s/Mpc     (Omega)mat = 0.2       (Omega)rad = 0.1

68  ENTER^
0.2  ENTER^
0.1  XEQ "AUL"   >>>>    T  = 10.35899203  Gy
RDN     P  = 0.000000000  Gy   ( no periodic )
RDN    R0 =  17.18654760  Gly
RDN    k  =  -1  ( hyperbolic space )

and    R05 = 0.000000000 Gly     = minimum radius of this Universe
R07 = 9.9999999 E99 Gly = maximum radius of this Universe = infinity

Example3:    H0 = 68 km/s/Mpc     (Omega)mat = 1.2       (Omega)rad = -2

68  ENTER^
1.2  ENTER^
-2   XEQ "AUL"   >>>>    T  =   5.732772984  Gy
RDN     P  = 0.000000000  Gy   ( no periodic )
RDN    R0 =  10.71769547  Gly
RDN    k  =  -1  ( hyperbolic space )

and    R05 = 8.276293007 Gly     = minimum radius of this Universe
R07 = 9.9999999 E99 Gly = maximum radius of this Universe = infinity

Example4:    H0 = 68 km/s/Mpc     (Omega)mat = 1.2       (Omega)rad = -0.2

68  ENTER^
1.2  ENTER^
-0.2  XEQ "AUL"   >>>>    T  = 10.65133139  Gy
RDN     P  = 0.000000000  Gy   ( no periodic )
RDN    R0 =  14.37929737  Gly
RDN    k  =  0  ( Euclidean space )

and    R05 = 2.396549562 Gly     = minimum radius of this Universe
R07 = 9.9999999 E99 Gly = maximum radius of this Universe = infinity

Example5:    H0 = 68 km/s/Mpc     (Omega)mat = 1       (Omega)rad = 0

68  ENTER^
1   ENTER^
0  XEQ "AUL"   >>>>    T  = 9.586198246  Gy
RDN     P  = 0.000000000  Gy   ( no periodic )
RDN    R0 =  14.37929737  Gly
RDN    k  =  0  ( Euclidean space )

and    R05 = 0.000000000 Gly     = minimum radius of this Universe
R07 = 9.9999999 E99 Gly = maximum radius of this Universe = infinity

Example6:    H0 = 68 km/s/Mpc     (Omega)mat = 0.4       (Omega)rad = 0.6

68   ENTER^
0.4   ENTER^
0.6  XEQ "AUL"   >>>>   T  = 7.759788005  Gy
RDN     P  = 0.000000000  Gy   ( no periodic )
RDN    R0 =  14.37929737  Gly
RDN    k  =  0  ( Euclidean space )

and    R05 = 0.000000000 Gly     = minimum radius of this Universe
R07 = 9.9999999 E99 Gly = maximum radius of this Universe = infinity

Example7:    H0 = 68 km/s/Mpc     (Omega)mat = 0       (Omega)rad = 1

68   ENTER^
0   ENTER^
1  XEQ "AUL"   >>>>   T  = 7.189648685  Gy
RDN     P  = 0.000000000  Gy   ( no periodic )
RDN    R0 =  14.37929737  Gly
RDN    k  =  0  ( Euclidean space )

and    R05 = 0.000000000 Gly     = minimum radius of this Universe
R07 = 9.9999999 E99 Gly = maximum radius of this Universe = infinity

Example8:    H0 = 68 km/s/Mpc     (Omega)mat = -0.4       (Omega)rad = 1.4

68   ENTER^
-0.4   ENTER^
1.4  XEQ "AUL"   >>>>   T  = 6.770532843  Gy
RDN     P  = 396.9889956  Gy
RDN    R0 =  14.37929737  Gly
RDN    k  =  0  ( Euclidean space )

and    R05 = 0.000000000 Gly     = minimum radius of this Universe
R07 = 50.32754080 Gly = maximum radius of this Universe

Example9:    H0 = 68 km/s/Mpc     (Omega)mat = -0.4       (Omega)rad = 1.3999

68     ENTER^
-0.4     ENTER^
1.3999  XEQ "AUL"   >>>>    T  = 6.770851649  Gy
RDN     P  = 397.7828431  Gy
RDN    R0 =  1437.929737  Gly
RDN    k  =  -1  ( hyperbolic space )

and    R05 =  0                   Gly  = minimum radius of this Universe
R07 = 5036.805349 Gly = maximum radius of this Universe.

Notes:

-In examples 1 to 8, the  results have a good precision.
-But in example 9, free42 gives   T = 6.770673328 & P = 397.7820792
-When  1 - (Omega)mat - (Omega)rad  is very small but # 0, this program does not give accurate results with an HP41 ( the outputs T & P may even be meaningless... )

b)  Cosmological Constant # 0

-This program employs Romberg integration.
-The precision is controlled by the display format ( lines 395...398 ).

Data Registers:    R00 = H0

R01 =
(Omega)mat            R04 = T          R06 = Rmin         R08 = R0      R10.....:  temp
R02 =
(Omega)lambda       R05 = P          R07 = Rmax         R09 = k

Flags:  F08-F09-F10-F24
Subroutines: /

-Lines 152-230-239-248 are three-byte GTOs

 01 LBL "AUM"  02 CF 08  03 CF 09  04 CF 10  05 R^  06 STO 00  07 SIGN  08 X<>Y  09 STO 03  10 STO 05  11 -  12 X<>Y  13 STO 02  14 -  15 X<>Y  16 STO 01  17 STO 04  18 -  19 STO 10  20 RCL 02  21 ST/ 04  22 ST/ 05  23 ST/ 10  24 CLX  25 STO 06  26 STO 07  27 LBL 01  28 RCL 06  29 X^2  30 STO 11  31 RCL 10  32 +  33 STO 13  34 RCL 07  35 ST+ X  36 STO 12  37 -  38 STO 08  39 RCL 04  40 RCL 06  41 RCL 08  42 *  43 -  44 STO 09  45 *  46 RCL 07  47 RCL 13  48 -  49 RCL 07  50 *  51 RCL 05  52 +  53 STO 14  54 RCL 06  55 ST+ X  56 STO 15  57 *  58 +  59 RCL 10  60 RCL 11  61 3  62 *  63 +  64 RCL 12  65 -  66 ST* 08  67 RCL 14            68 *  69 RCL 06  70 RCL 09  71 *  72 RCL 12 73 *  74 -  75 RCL 07  76 RCL 15  77 X^2  78 *  79 RCL 08  80 +  81 ST/ Z  82 /  83 ST+ 07  84 X<>Y  85 ST+ 06  86 R-P  87 VIEW X  88  E-8  89 XY 123 STO 07 124 LBL 03 125 RCL 09 126 RCL 08 127 2 128 / 129 CHS 130 ENTER 131 X^2 132 RCL 09 133 - 134 X<0? 135 SF 10 136 X<0? 137 GTO 07 138 SQRT 139 RCL Y 140 SIGN 141 * 142 + 143 X#0? 144 ST/ Y 145 STO 08 146 X<>Y 147 STO 09 148 LBL 07 149 FS? 09 150 FC? 10 151 FS? 30 152 GTO 13 153 FC? 09 154 GTO 03 155 RCL 08 156 X<> 06 157 STO 08 158 RCL 09 159 X<> 07 160 STO 09 161 LBL 03 162 FC?C 10 163 FS? 09 164 GTO 04 165 RCL 06 166 RCL 07 167 X>Y? 168 X<>Y 169 RCL 08 170 X>Y? 171 X<>Y 172 RCL 09 173 X>Y? 174 X<>Y 175 STO 06 176 RDN 177 XY 179 RDN 180 X>Y? 181 X<>Y 182 STO 07 183 X<> T 184 XY 186 STO 09 187 X<>Y 188 STO 08 189 1 190 X>Y? 191 GTO 03 192 RCL 07         193 X>Y? 194 GTO 11 195 X<> 06 196 X<> 08 197 STO 07 198 GTO 11 199 LBL 03 200 RCL 08 201 X<> 06 202 STO 08 203 RCL 09 204 X<> 07 205 STO 09 206 LBL 11 207 RCL 08 208 RCL 09 209 + 210 CHS 211 X<> 08 212 ST* 09 213 GTO 07 214 LBL 04 215 RCL 06 216 RCL 07 217 X>Y? 218 X<>Y 219 STO 06 220 X<>Y 221 STO 07 222 CF 09 223 LBL 07 224 RCL 06 225 1 226 XY 284 RCL 06 285 - 286 * 287 SQRT 288 / 289 RTN 290 LBL 09 291 ENTER 292 ENTER 293 X^2 294 RCL 10 295 + 296 * 297 RCL 04 298 + 299 * 300 RCL 05 301 + 302 SQRT 303 / 304 RTN 305 LBL 11 306 X^2 307 RCL 06 308 X<>Y 309 - 310 ENTER 311 STO Z 312 RCL 08 313 + 314 * 315 RCL 09 316 + 317 RCL 07 318 R^ 319 - 320 * 321 SQRT 322 / 323 CHS 324 RTN 325 LBL 10 326 STO 11 327 X<>Y 328 STO 12 329 STO 13 330 XEQ 00 331 STO 15 332 RCL 11 333 ST- 13 334 XEQ 00 335 RCL 15 336 + 337 2 338 / 339 STO 15 340 RCL 13 341 * 342 STO 20 343 20 344 STO 19 345 SIGN 346 STO 17 347 LBL 02 348 RCL 17 349 STO 16 350 RCL 11 351 RCL 13 352 2 353 / 354 + 355 LBL 05 356 STO 14         357 XEQ 00 358 ST+ 15 359 RCL 13 360 RCL 14 361 + 362 DSE 16 363 GTO 05 364 2 365 ST/ 13 366 ST* 17 367 X^2 368 STO 18 369 RCL 19 370 INT 371  E3 372 / 373 20 374 + 375 STO 19 376 RCL 13 377 RCL 15 378 * 379 LBL 06 380 ENTER 381 ENTER 382 X<> IND 19 383 STO 14 384 - 385 RCL 18 386 4 387 ST* 18 388 SIGN 389 - 390 / 391 + 392 ISG 19 393 GTO 06 394 STO IND 19 395 VIEW X 396 RND 397 RCL 14 398 RND 399 X#Y? 400 GTO 02 401 RCL IND 19 402 RTN 403 LBL 04 404 - 405 SQRT 406 RCL 07 407 CHS 408 X<0? 409 CLX 410 SF 10 411 SQRT 412 XEQ 10 413 ST+ X 414 STO 04 415 RCL 07 416 STO 06 417 CLX 418 STO 07 419 GTO 03 420 LBL 13 421 CLX 422 SIGN 423 LASTX 424 XEQ 10 425 STO 04 426 CLX 427 STO 07         428 GTO 03 429 LBL 12 430 XEQ 10 431 STO 04 432 CLST 433 RCL 12 434 XEQ 10 435 RCL 06 436 STO 07 437 CLX 438 STO 06 439 X<>Y 440 GTO 12 441 LBL 14 442 RCL 06 443 RCL 06 444 RCL 07 445 - 446 / 447 X<0? 448 CLX 449 SQRT 450 ASIN 451 XEQ 10 452 STO 04 453 SIGN 454 ASIN 455 RCL 12 456 XEQ 10 457 LBL 12 458 RCL 04 459 ST+ 04 460 + 461 4 462 * 463 LBL 03 464 STO 05 465 RCL 06 466 X<0? 467 CLX 468 STO 06 469 DEG 470 90 471 TAN 472 RCL 07 473 X<=0? 474 X<>Y 475 STO 07 476 SF 24 477 RCL 02 478 ABS 479 SQRT 480 ST/ 04 481 ST/ 05 482 RCL 02 483 RCL 10 484 * 485 STO 14 486 CHS 487 X#0? 488 SIGN 489 STO 09 490 977.7922214 491 RCL 00 492 / 493 ST* 04 494 ST* 05 495 RCL 14         496 ABS 497 SQRT 498 X=0? 499 SIGN 500 / 501 ST* 06 502 ST* 07 503 STO 08 504 RCL 05 505 RCL 04 506 CLD 507 END

( 649 bytes / SIZE 020+??? )

 STACK INPUTS OUTPUTS T H0 > 0 k Z (Omega) mat R0 Y (Omega)lambda P X (Omega)rad T

T = Age of the Universe in Gigayears  ,       k = +1 , 0 , -1  =  spherical, euclidean , hyperbolic Universes
P = Period of the Universe in Gigayears   ,   R0 = current scale factor of the Universe in Gigalightyears

Example1:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001

FIX 7

71    ENTER^
-0.1   ENTER^
-0.003  ENTER^
-0.001  XEQ "AUM"   >>>>   T  = 14.61256011                            ---Execution time = 2m52s---
RDN   P =  796.1455562
RDN   R0 = 13.10701187
RDN    k = -1 ( hyperbolic space )

-We also have in R06 & R07

Rmin = 1.306380222
Rmax = 250.8400402

Example2:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049   ( our Universe ? )

FIX 7

71      ENTER^
0.044    ENTER^
0.521    ENTER^
0.000049  XEQ "AUM"   >>>>   T  = 15.60221461                            ---Execution time = 2m37s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 20.88180628
RDN    k = -1

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example3:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,   (Omega) rad  = 0.000049   ( our Universe ? )

FIX 7

71       ENTER^
0.271     ENTER^
0.732     ENTER^
0.000049   XEQ "AUM"   >>>>   T  = 13.66788682                            ---Execution time = 2m31s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 249.4075046
RDN    k = +1   ( spherical space )

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example4:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.6 ,  (Omega)lambda  = 0.4 ,   (Omega) rad  = 0.1

FIX 7

71     ENTER^
-0.6    ENTER^
0.4    ENTER^
0.1    XEQ "AUM"   >>>>   T  = 17.27498133                             ---Execution time = 4m13s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 13.13082118
RDN    k = -1   ( hyperbolic space )

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example5:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.01 ,  (Omega)lambda  = -0.4 ,   (Omega) rad  = 0.01

FIX 7

71    ENTER^
0.01   ENTER^
-0.4   ENTER^
0.01   XEQ "AUM"   >>>>   T  = 11.30457454                            ---Execution time = 7m18s---
RDN   P =  66.21313843
RDN   R0 = 11.72326781
RDN    k = -1 ( hyperbolic space )

Rmin = 0
Rmax = 21.83998798

Example6:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.4 ,  (Omega)lambda  = +0.001 ,   (Omega) rad  = 1.4

FIX 7

71    ENTER^
-0.4   ENTER^
0.001  ENTER^
1.4    XEQ "AUM"   >>>>   T  = 6.484937796                            ---Execution time = 2m32s---
RDN   P =  900.0920456
RDN   R0 = 435.5000702
RDN    k = +1 ( spherical space )

Rmin = 0
Rmax = 1889.773264

Notes:

-We use Bairstow method ( LBL 01 ) to factorize the quartic polynomial before finding its roots.
-Romberg method = LBL 10
-FIX 9 would often produce infinite loops

-This program also works if (Omega)mat = 0 , but in this case, "AUR" listed in paragraph 2°) is much faster !
-It doesn't work if (Omega)lambda  = 1 and  (Omega)mat = 0 = (Omega) rad but in this case,  T = +infinity ,  Rmin = 0  and   Rmax = +infinity  ( de Sitter model )
-We could add a few lines to take this case into account...

-We can save several bytes if we use Gauss-Legendre 2-point formula.

-The first estimations are computed with 15 subintervals for T and 55 subintervals for P ( lines 414 & 415 )
-After these results, you can use other values  N1 for T & N2  for P, for instance, with  30 & 100

30   ENTER^
100  XEQ 10

Data Registers:    R00 = H0

R01 =
(Omega)mat            R21 = T          R23 = Rmin         R25 = R0
R02 =
(Omega)lambda       R22 = P          R24 = Rmax         R26 = k

Flags:  F08-F09-F10-F24
Subroutines: /

-Lines 152-230-239-250 are three-byte GTOs

 01 LBL "AUM"  02 CF 08  03 CF 09  04 CF 10  05 R^  06 STO 00  07 SIGN  08 X<>Y  09 STO 03  10 STO 05  11 -  12 X<>Y  13 STO 02  14 -  15 X<>Y  16 STO 01  17 STO 04  18 -  19 STO 10  20 RCL 02  21 ST/ 04  22 ST/ 05  23 ST/ 10  24 CLX  25 STO 06  26 STO 07  27 LBL 01  28 RCL 06  29 X^2  30 STO 11  31 RCL 10  32 +  33 STO 13  34 RCL 07  35 ST+ X  36 STO 12  37 -  38 STO 08  39 RCL 04  40 RCL 06  41 RCL 08  42 *  43 -  44 STO 09  45 *  46 RCL 07  47 RCL 13  48 -  49 RCL 07  50 *  51 RCL 05  52 +  53 STO 14  54 RCL 06  55 ST+ X  56 STO 15  57 *  58 +  59 RCL 10  60 RCL 11  61 3  62 *  63 +  64 RCL 12            65 -  66 ST* 08  67 RCL 14  68 * 69 RCL 06  70 RCL 09  71 *  72 RCL 12  73 *  74 -  75 RCL 07  76 RCL 15  77 X^2  78 *  79 RCL 08  80 +  81 ST/ Z  82 /  83 ST+ 07  84 X<>Y  85 ST+ 06  86 R-P  87 VIEW X  88  E-8  89 XY 123 STO 07 124 LBL 03 125 RCL 09 126 RCL 08 127 2 128 / 129 CHS 130 ENTER 131 X^2 132 RCL 09         133 - 134 X<0? 135 SF 10 136 X<0? 137 GTO 07 138 SQRT 139 RCL Y 140 SIGN 141 * 142 + 143 X#0? 144 ST/ Y 145 STO 08 146 X<>Y 147 STO 09 148 LBL 07 149 FS? 09 150 FC? 10 151 FS? 30 152 GTO 13 153 FC? 09 154 GTO 03 155 RCL 08 156 X<> 06 157 STO 08 158 RCL 09 159 X<> 07 160 STO 09 161 LBL 03 162 FC?C 10 163 FS? 09 164 GTO 04 165 RCL 06 166 RCL 07 167 X>Y? 168 X<>Y 169 RCL 08 170 X>Y? 171 X<>Y 172 RCL 09 173 X>Y? 174 X<>Y 175 STO 06 176 RDN 177 XY 179 RDN 180 X>Y? 181 X<>Y 182 STO 07 183 X<> T 184 XY 186 STO 09 187 X<>Y 188 STO 08 189 1 190 X>Y? 191 GTO 03 192 RCL 07 193 X>Y? 194 GTO 11 195 X<> 06 196 X<> 08 197 STO 07         198 GTO 11 199 LBL 03 200 RCL 08 201 X<> 06 202 STO 08 203 RCL 09 204 X<> 07 205 STO 09 206 LBL 11 207 RCL 08 208 RCL 09 209 + 210 CHS 211 X<> 08 212 ST* 09 213 GTO 07 214 LBL 04 215 RCL 06 216 RCL 07 217 X>Y? 218 X<>Y 219 STO 06 220 X<>Y 221 STO 07 222 CF 09 223 LBL 07 224 RCL 06 225 1 226 XY 255 STO 11 256 ST+ 11 257 / 258 STO Y 259 3 260 SQRT 261 / 262 STO 13 263 - 264 STO 14         265 2 266 / 267 ST- 12 268 CLX 269 STO 15 270 LBL 05 271 RCL 12 272 RCL 13 273 X<> 14 274 STO 13 275 + 276 STO 12 277 XEQ 00 278 ST+ 15 279 DSE 11 280 GTO 05 281 RCL 13 282 RCL 15 283 * 284 3 285 SQRT 286 * 287 RTN 288 LBL 00 289 FS? 08 290 GTO 11 291 FS? 09 292 GTO 09 293 FS? 10 294 GTO 08 295 SIN 296 X^2 297 RCL 05 298 * 299 RCL 06 300 + 301 ENTER 302 STO Z 303 RCL 08 304 + 305 * 306 RCL 09 307 + 308 SQRT 309 / 310 RTN 311 LBL 08 312 X^2 313 RCL 07 314 + 315 ENTER 316 STO Z 317 RCL 08 318 + 319 * 320 RCL 09 321 + 322 X<>Y 323 RCL 06 324 - 325 * 326 SQRT 327 / 328 RTN 329 LBL 09 330 ENTER 331 ENTER 332 X^2 333 RCL 10         334 + 335 * 336 RCL 04 337 + 338 * 339 RCL 05 340 + 341 SQRT 342 ST+ X 343 / 344 RTN 345 LBL 11 346 X^2 347 RCL 06 348 X<>Y 349 - 350 ENTER 351 STO Z 352 RCL 08 353 + 354 * 355 RCL 09 356 + 357 RCL 07 358 R^ 359 - 360 * 361 SQRT 362 / 363 CHS 364 RTN 365 LBL 04 366 - 367 SQRT 368 STO 17 369 RCL 07 370 STO 23 371 CHS 372 X<0? 373 CLX 374 SF 10 375 SQRT 376 STO 16 377 CLX 378 STO 24 379 GTO 04 380 LBL 13 381 CLX 382 STO 16 383 STO 24 384 SIGN 385 STO 17 386 RCL 06 387 STO 23 388 GTO 04 389 LBL 12 390 CLX 391 STO 18 392 STO 23 393 RCL 06 394 STO 24 395 GTO 04 396 LBL 14 397 STO 17 398 RCL 06 399 RCL 06 400 STO 23         401 RCL 07 402 STO 24 403 - 404 / 405 X<0? 406 CLX 407 SQRT 408 ASIN 409 STO 16 410 SIGN 411 ASIN 412 STO 18 413 LBL 04 414 15 415 55 416 LBL 10 417 RAD 418 STO 20 419 X<>Y 420 STO 19 421 RCL 17 422 RCL 16 423 XEQ 06 424 STO 21 425 CLX 426 FC? 09 427 FS? 10 428 GTO 00 429 RCL 20 430 RCL 18 431 RCL 16 432 XEQ 06 433 ST+ X 434 LBL 00 435 STO 22 436 RCL 23 437 X<0? 438 CLX 439 STO 23 440 DEG 441 90 442 TAN 443 RCL 24 444 X<=0? 445 X<>Y 446 STO 24 447 SF 24 448 RCL 02 449 ABS 450 SQRT 451 ST/ 21 452 ST/ 22 453 RCL 02 454 RCL 10 455 * 456 STO 14 457 CHS 458 X#0? 459 SIGN 460 STO 26 461 977.7922214 462 RCL 00 463 / 464 ST* 21 465 ST* 22 466 RCL 14         467 ABS 468 SQRT 469 X=0? 470 SIGN 471 / 472 ST* 23 473 ST* 24 474 STO 25 475 RCL 22 476 RCL 21 477 CF 24 478 END

( 631 bytes / SIZE 027 )

 STACK INPUTS OUTPUTS T H0 > 0 k Z (Omega) mat R0 Y (Omega)lambda P X (Omega)rad T

T = Age of the Universe in Gigayears  ,       k = +1 , 0 , -1  =  spherical, euclidean , hyperbolic Universes
P = Period of the Universe in Gigayears   ,   R0 = current scale factor of the Universe in Gigalightyears

Example1:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001

71    ENTER^
-0.1   ENTER^
-0.003  ENTER^
-0.001  XEQ "AUM"   >>>>   T  = 14.61256012                            ---Execution time = 3m33s---
RDN   P =  796.1455560
RDN   R0 = 13.10701187
RDN    k = -1 ( hyperbolic space )

-We also have in R23 & R24

Rmin = 1.306380222
Rmax = 250.8400402

Example2:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049   ( our Universe ? )

71      ENTER^
0.044    ENTER^
0.521    ENTER^
0.000049  XEQ "AUM"   >>>>   T  = 15.60221156                            ---Execution time = 48s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 20.88180628
RDN    k = -1

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

-With another number of subintervals to compute T - for instance 30 -

30   ENTER^  ( any other number in X-register since P = 0 )   XEQ 10  >>>>   T= 15.60221435

-With N = 60 , it yields:   T = 15.60221460

... and so on ...

Note:

-Of course, we can use Carlson elliptic integrals to obtain very accurate results:

Data Registers:    R00 = H0

R01 =
(Omega)mat            R04 = T          R06 = Rmin         R08 = R0      R10.....:  temp
R02 =
(Omega)lambda       R05 = P          R07 = Rmax         R09 = k

Flags:  F08-F09-F10-F24
Subroutines:  ( M-Code routines )  RF  RFZ  RJ  RJZ  ( cf "Carlson Elliptic Integrals for the HP41" )

-Lines 156-334 are three-byte GTOs

 01 LBL "AUM"  02 CF 08  03 CF 09  04 CF 10  05 R^  06 STO 00  07 SIGN  08 X<>Y  09 STO 03  10 STO 05  11 -  12 X<>Y  13 STO 02  14 -  15 X<>Y  16 STO 01  17 STO 04  18 -  19 STO 10  20 STO 20  21 RCL 02  22 ST/ 04  23 ST/ 05  24 ST/ 10  25 CLX  26 STO 06  27 STO 07  28 LBL 01  29 RCL 06  30 X^2  31 STO 11  32 RCL 10  33 +  34 STO 13  35 RCL 07  36 ST+ X  37 STO 12  38 -  39 STO 08  40 RCL 04  41 RCL 06  42 RCL 08  43 *  44 -  45 STO 09  46 *  47 RCL 07  48 RCL 13  49 -  50 RCL 07  51 *  52 RCL 05  53 +  54 STO 14  55 RCL 06  56 ST+ X  57 STO 15  58 *  59 +  60 RCL 10  61 RCL 11  62 3  63 *  64 +  65 RCL 12  66 -  67 ST* 08  68 RCL 14  69 *  70 RCL 06  71 RCL 09  72 *  73 RCL 12  74 *  75 -  76 RCL 07  77 RCL 15  78 X^2  79 *  80 RCL 08  81 +  82 ST/ Z  83 /  84 ST+ 07  85 X<>Y  86 ST+ 06  87 R-P  88 VIEW X  89  E-8  90 XY 126 STO 07 127 LBL 02 128 RCL 09 129 RCL 08 130 2 131 / 132 CHS 133 STO 13 134 ENTER 135 X^2 136 RCL 09 137 - 138 STO 14 139 X<0? 140 SF 10 141 X<0? 142 GTO 02 143 SQRT 144 RCL Y 145 SIGN 146 * 147 + 148 X#0? 149 ST/ Y 150 STO 08 151 X<>Y 152 STO 09 153 LBL 02 154 FS? 09 155 FC? 10 156 GTO 03 157 1 158 RCL 07 159 + 160 RCL 06 161 + 162 SQRT 163 STO 11 164 SIGN 165 RCL 06 166 - 167 RCL 09 168 + 169 SQRT 170 STO 12 171 RCL 09 172 SQRT 173 + 174 X^2 175 1 176 - 177 SQRT 178 RCL 11 179 RCL 07 180 SQRT 181 + 182 X^2 183 1 184 - 185 SQRT 186 * 187 X^2 188 STO 15 189 RCL 07 190 RCL 09 191 + 192 ST+ X 193 RCL 06 194 X^2 195 STO 14 196 + 197 STO 16 198 X^2 199 RCL 07 200 4 201 * 202 RCL 14 203 - 204 STO 08 205 RCL 09 206 4 207 * 208 RCL 14 209 - 210 * 211 - 212 SQRT 213 ENTER 214 CHS 215 RCL 15 216 RCL 16 217 + 218 STO 18 219 ST+ Z 220 + 221 STO 14           222 X<>Y 223 STO 13 224 RCL 15 225 RCL 08 226 + 227 STO 17 228 X<> Z 229 RCL 15 230 RJ 231 3 232 / 233 X<> 18 234 2 235 / 236 RCL 11 237 RCL 09 238 SQRT 239 * 240 RCL 12 241 RCL 07 242 SQRT 243 * 244 + 245 STO 19 246 X^2 247 - 248 STO 16 249 RCL 17 250 * 251 RCL 19 252 / 253 RCL 08 254 ST+ X 255 RCL 19 256 * 257 + 258 X^2 259 RCL 16 260 RCL 19 261 / 262 X^2 263 RCL 08 264 + 265 RCL 17 266 X^2 267 * 268 ENTER 269 RF 270 2 271 / 272 RCL 18 273 + 274 RCL 08 275 * 276 RCL 06 277 ST+ X 278 * 279 STO 18 280 RCL 16 281 RCL 11 282 ST/ 12 283 RCL 07 284 SQRT 285 * 286 STO 11 287 / 288 2 289 + 290 X^2 291 RCL 16 292 X^2 293 RCL 19 294 X^2 295 RCL 08 296 * 297 + 298 RCL 11 299 X^2 300 / 301 ENTER 302 RF 303 RCL 06 304 2 305 + 306 RCL 12 307 * 308 RCL 09 309 RCL 07 310 / 311 SQRT 312 RCL 06 313 * 314 + 315 * 316 2 317 / 318 ST+ 18 319 RCL 15 320 RCL 14 321 RCL 13 322 RF 323 RCL 06 324 * 325 CHS 326 RCL 18         327 + 328 ST+ X 329 STO 04 330 CLX 331 STO 05 332 STO 06 333 STO 07 334 GTO 07 335 LBL 05 336 STO 11 337 RCL 06 338 - 339 RCL 11 340 RCL 07 341 - 342 FS? 09 343 CHS 344 STO 13 345 * 346 X<>Y 347 STO 12 348 RCL 08 349 - 350 FS? 10 351 GTO 00 352 * 353 RCL 12 354 RCL 09 355 - 356 * 357 SQRT 358 RCL 12 359 RCL 06 360 - 361 RCL 12 362 RCL 07 363 - 364 FS? 09 365 CHS 366 STO 14 367 * 368 RCL 11 369 RCL 08 370 - 371 * 372 RCL 11 373 RCL 09 374 - 375 * 376 SQRT 377 + 378 X^2 379 STO 15 380 STO 18 381 RCL 11 382 RCL 06 383 - 384 RCL 11 385 RCL 08 386 - 387 * 388 RCL 14 389 * 390 RCL 12 391 RCL 09 392 - 393 * 394 SQRT 395 RCL 12 396 RCL 06 397 - 398 RCL 12 399 RCL 08 400 - 401 * 402 RCL 13 403 * 404 RCL 11 405 RCL 09 406 - 407 * 408 SQRT 409 + 410 X^2 411 STO 16 412 RCL 11 413 RCL 06 414 - 415 RCL 11 416 RCL 09 417 - 418 * 419 RCL 14 420 * 421 RCL 12 422 RCL 08 423 - 424 * 425 SQRT 426 RCL 12 427 RCL 06 428 - 429 RCL 12 430 RCL 09 431 - 432 * 433 RCL 13 434 * 435 RCL 11         436 RCL 08 437 - 438 * 439 SQRT 440 + 441 X^2 442 STO 17 443 GTO 02 444 LBL 00 445 X^2 446 RCL 09 447 X^2 448 + 449 * 450 SQRT 451 RCL 12 452 RCL 06 453 - 454 RCL 12 455 RCL 07 456 - 457 FS? 09 458 CHS 459 STO 14 460 * 461 RCL 11 462 RCL 08 463 - 464 X^2 465 RCL 09 466 X^2 467 + 468 * 469 SQRT 470 + 471 X^2 472 STO 15 473 STO 18 474 RCL 11 475 RCL 12 476 - 477 RCL 09 478 * 479 STO 05 480 RCL 11 481 RCL 08 482 - 483 RCL 12 484 LASTX 485 - 486 * 487 RCL 09 488 X^2 489 + 490 STO 19 491 RCL 11 492 RCL 06 493 - 494 RCL 14 495 * 496 ST* Z 497 * 498 R-P 499 X<>Y 500 2 501 / 502 X<>Y 503 SQRT 504 P-R 505 STO 16 506 X<>Y 507 STO 17 508 RCL 05 509 CHS 510 RCL 19 511 RCL 12 512 RCL 06 513 - 514 RCL 13 515 * 516 ST* Z 517 * 518 R-P 519 X<>Y 520 2 521 / 522 X<>Y 523 SQRT 524 P-R 525 X<>Y 526 RCL 17 527 + 528 X<>Y 529 RCL 16 530 + 531 R-P 532 X<>Y 533 ST+ X 534 X<>Y 535 X^2 536 P-R 537 STO 16 538 X<>Y 539 STO 17 540 LBL 02 541 RCL 11 542 RCL 12 543 - 544 X^2 545 ST/ 15 546 ST/ 16 547 ST/ 17 548 ST/ 18 549 RCL 09 550 RCL 06         551 - 552 RCL 08 553 LASTX 554 - 555 ST* Y 556 X^2 557 RCL 09 558 X^2 559 + 560 FC? 10 561 X<>Y 562 STO 19 563 FC? 09 564 CHS 565 ST+ 18 566 RCL 18 567 RCL 17 568 RCL 16 569 RCL 15 570 FS? 10 571 RJZ 572 FC? 10 573 RJ 574 RCL 06 575 RCL 07 576 - 577 FC? 08 578 FS? 09 579 CHS 580 * 581 3 582 / 583 ST* 19 584 RCL 18 585 RCL 11 586 RCL 06 587 - 588 RCL 12 589 LASTX 590 - 591 * 592 X=0? 593 GTO 00 594 / 595 ENTER 596 STO Z 597 1 598 FS? 09 599 CHS 600 + 601 XEQ "RF" 602 LBL 00 603 ST+ 19 604 RCL 17 605 RCL 16 606 RCL 15 607 FS? 10 608 RFZ 609 FC? 10 610 RF 611 RCL 06 612 * 613 RCL 19 614 FS? 08 615 CHS 616 + 617 ST+ X 618 RTN 619 LBL 03 620 FC? 09 621 GTO 03 622 RCL 08 623 STO 06 624 RCL 09 625 STO 07 626 RCL 11 627 STO 08 628 RCL 12 629 CHS 630 SQRT 631 STO 09 632 LBL 03 633 FC? 10 634 GTO 03 635 RCL 13 636 STO 08 637 RCL 14 638 CHS 639 SQRT 640 STO 09 641 LBL 03 642 FC? 10 643 FS? 09 644 GTO 04 645 RCL 06 646 RCL 07 647 X>Y? 648 X<>Y 649 RCL 08 650 X>Y? 651 X<>Y 652 RCL 09 653 X>Y? 654 X<>Y 655 STO 06 656 RDN 657 XY 659 RDN 660 X>Y? 661 X<>Y 662 STO 07 663 X<> T 664 XY 666 STO 09         667 X<>Y 668 STO 08 669 1 670 X>Y? 671 GTO 02 672 RCL 07 673 X>Y? 674 GTO 03 675 X<> 06 676 X<> 08 677 STO 07 678 GTO 03 679 LBL 02 680 RCL 08 681 X<> 06 682 STO 08 683 RCL 09 684 X<> 07 685 STO 09 686 GTO 03 687 LBL 04 688 RCL 06 689 RCL 07 690 X>Y? 691 X<>Y 692 STO 06 693 X<>Y 694 STO 07 695 LBL 03 696 FS?C 09 697 SF 10 698 RCL 06 699 1 700 XY? 704 GTO 06 705 CF 09 706 RCL 06 707 X<> 07 708 STO 06 709 X<0? 710 0 711 ENTER 712 SIGN 713 XEQ 05 714 STO 04 715 CLX 716 STO 05 717 STO 07 718 GTO 07 719 LBL 02 720 CF 09 721 SF 08 722 CLST 723 SIGN 724 XEQ 05 725 STO 04 726 CLST 727 RCL 06 728 XEQ 05 729 ST+ X 730 STO 05 731 CLX 732 X<> 06 733 STO 07 734 GTO 07 735 LBL 06 736 SF 09 737 RCL 06 738 X<0? 739 0 740 ENTER 741 SIGN 742 XEQ 05 743 STO 04 744 RCL 06 745 X<0? 746 0 747 RCL 07 748 XEQ 05 749 ST+ X 750 STO 05 751 LBL 07 752 RCL 06 753 X<0? 754 CLX 755 STO 06 756 DEG 757 90 758 TAN 759 RCL 07 760 X<=0? 761 X<>Y 762 STO 07 763 SF 24 764 RCL 02 765 ABS 766 SQRT 767 ST/ 04 768 ST/ 05 769 RCL 20 770 CHS 771 X#0? 772 SIGN 773 STO 09 774 977.7922214 775 RCL 00         776 / 777 ST* 04 778 ST* 05 779 RCL 20 780 ABS 781 SQRT 782 X=0? 783 SIGN 784 / 785 ST* 06 786 ST* 07 787 STO 08 788 RCL 05 789 RCL 04 790 CF 24 791 END

( 977 bytes / SIZE 021 )

 STACK INPUTS OUTPUTS T H0 > 0 k Z (Omega) mat R0 Y (Omega)lambda P X (Omega)rad T

T = Age of the Universe in Gigayears  ,       k = +1 , 0 , -1  =  spherical, euclidean , hyperbolic Universes
P = Period of the Universe in Gigayears   ,   R0 = current scale factor of the Universe in Gigalightyears

Example1:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001

71    ENTER^
-0.1   ENTER^
-0.003  ENTER^
-0.001  XEQ "AUM"   >>>>   T  = 14.61256009                           ---Execution time = 59s---
RDN   P =  796.1455565
RDN   R0 = 13.10701187
RDN    k = -1 ( hyperbolic space )

-We also have in R06 & R07

Rmin = 1.306380222
Rmax = 250.8400402

Example2:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049   ( our Universe ? )

71      ENTER^
0.044    ENTER^
0.521    ENTER^
0.000049  XEQ "AUM"   >>>>   T  = 15.60221452                            ---Execution time = 48s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 20.88180628
RDN    k = -1

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example3:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,   (Omega) rad  = 0.000049   ( our Universe ? )

71       ENTER^
0.271     ENTER^
0.732     ENTER^
0.000049   XEQ "AUM"   >>>>   T  = 13.66788681                            ---Execution time = 102s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 249.4075046
RDN    k = +1   ( spherical space )

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example4:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.6 ,  (Omega)lambda  = 0.4 ,   (Omega) rad  = 0.1

71     ENTER^
-0.6    ENTER^
0.4    ENTER^
0.1    XEQ "AUM"   >>>>   T  = 17.27498133                             ---Execution time = 44s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 13.13082118
RDN    k = -1   ( hyperbolic space )

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example5:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.01 ,  (Omega)lambda  = -0.4 ,   (Omega) rad  = 0.01

71    ENTER^
0.01   ENTER^
-0.4   ENTER^
0.01   XEQ "AUM"   >>>>   T  = 11.30457447                            ---Execution time = 90s---
RDN   P =  66.21313811
RDN   R0 = 11.72326781
RDN    k = -1 ( hyperbolic space )

Rmin = 0
Rmax = 21.83998798

Example6:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.4 ,  (Omega)lambda  = +0.001 ,   (Omega) rad  = 1.4

71    ENTER^
-0.4   ENTER^
0.001  ENTER^
1.4    XEQ "AUM"   >>>>   T  = 6.484937782                           ---Execution time = 114s---
RDN   P =  900.0920456
RDN   R0 = 435.5000702
RDN    k = +1 ( spherical space )

Rmin = 0
Rmax = 1889.773264

Notes:

-Lines 157 to 330 employ the formula given in reference [3]
-Lines 335 to 618 use the formula given in reference [1]         ( for elliptic integrals [-1,-1,-1;-1,+2] )

-But the integral:     §ar x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ]      may be re-written:

§0+infinity  sqrt(r-a)  (a.y+r)  /  (y+1) / sqrt [ ((b-a).y + (b-r)) ((a+c).y + (r+c)) ((a+d).y + (r+d)) ]  dy    by the change of variable:    x = (a.y+r) / (y+1)

-This reduces the number of bytes ( 758 instead of 976 ) , but the execution times will be a little increased... ( except in examples 1 & 4 )

Data Registers:    R00 = H0

R01 =
(Omega)mat            R04 = T          R06 = Rmin         R08 = R0      R10.....:  temp
R02 =
(Omega)lambda       R05 = P          R07 = Rmax         R09 = k

Flags:  F09-F10-F24
Subroutines:  ( M-Code routines )  RF  RFZ  RJ  RJZ  ( cf "Carlson Elliptic Integrals for the HP41" )

-Lines 155-333 are three-byte GTOs

 01 LBL "AUM"  02 CF 09  03 CF 10  04 R^  05 STO 00  06 SIGN  07 X<>Y  08 STO 03  09 STO 05  10 -  11 X<>Y  12 STO 02  13 -  14 X<>Y  15 STO 01  16 STO 04  17 -  18 STO 10  19 STO 20  20 RCL 02  21 ST/ 04  22 ST/ 05  23 ST/ 10  24 CLX  25 STO 06  26 STO 07  27 LBL 01  28 RCL 06  29 X^2  30 STO 11  31 RCL 10  32 +  33 STO 13  34 RCL 07  35 ST+ X  36 STO 12  37 -  38 STO 08  39 RCL 04  40 RCL 06  41 RCL 08  42 *  43 -  44 STO 09  45 *  46 RCL 07  47 RCL 13  48 -  49 RCL 07  50 *  51 RCL 05  52 +  53 STO 14  54 RCL 06  55 ST+ X  56 STO 15  57 *  58 +  59 RCL 10  60 RCL 11  61 3  62 *  63 +  64 RCL 12  65 -  66 ST* 08  67 RCL 14  68 *  69 RCL 06  70 RCL 09            71 *  72 RCL 12  73 *  74 -  75 RCL 07  76 RCL 15  77 X^2  78 *  79 RCL 08  80 +  81 ST/ Z  82 /  83 ST+ 07  84 X<>Y  85 ST+ 06  86 R-P 87 VIEW X  88  E-8  89 XY 125 STO 07 126 LBL 02 127 RCL 09 128 RCL 08 129 2 130 / 131 CHS 132 STO 13 133 ENTER 134 X^2 135 RCL 09 136 - 137 STO 14 138 X<0? 139 SF 10 140 X<0? 141 GTO 02 142 SQRT 143 RCL Y 144 SIGN 145 * 146 + 147 X#0? 148 ST/ Y 149 STO 08 150 X<>Y 151 STO 09 152 LBL 02 153 FS? 09 154 FC? 10 155 GTO 03 156 1 157 RCL 07 158 + 159 RCL 06 160 + 161 SQRT 162 STO 11 163 SIGN 164 RCL 06 165 - 166 RCL 09 167 + 168 SQRT 169 STO 12 170 RCL 09         171 SQRT 172 + 173 X^2 174 1 175 - 176 SQRT 177 RCL 11 178 RCL 07 179 SQRT 180 + 181 X^2 182 1 183 - 184 SQRT 185 * 186 X^2 187 STO 15 188 RCL 07 189 RCL 09 190 + 191 ST+ X 192 RCL 06 193 X^2 194 STO 14 195 + 196 STO 16 197 X^2 198 RCL 07 199 4 200 * 201 RCL 14 202 - 203 STO 08 204 RCL 09 205 4 206 * 207 RCL 14 208 - 209 * 210 - 211 SQRT 212 ENTER 213 CHS 214 RCL 15 215 RCL 16 216 + 217 STO 18 218 ST+ Z 219 + 220 STO 14 221 X<>Y 222 STO 13 223 RCL 15 224 RCL 08 225 + 226 STO 17 227 X<> Z 228 RCL 15 229 RJ 230 3 231 / 232 X<> 18 233 2 234 / 235 RCL 11 236 RCL 09 237 SQRT 238 * 239 RCL 12 240 RCL 07 241 SQRT 242 * 243 + 244 STO 19 245 X^2 246 - 247 STO 16 248 RCL 17 249 * 250 RCL 19 251 / 252 RCL 08 253 ST+ X 254 RCL 19 255 * 256 + 257 X^2 258 RCL 16 259 RCL 19 260 / 261 X^2 262 RCL 08 263 + 264 RCL 17 265 X^2 266 * 267 ENTER 268 RF 269 2 270 / 271 RCL 18 272 + 273 RCL 08 274 * 275 RCL 06 276 ST+ X 277 * 278 STO 18 279 RCL 16 280 RCL 11 281 ST/ 12 282 RCL 07 283 SQRT 284 * 285 STO 11 286 / 287 2 288 + 289 X^2 290 RCL 16 291 X^2 292 RCL 19 293 X^2 294 RCL 08 295 * 296 + 297 RCL 11 298 X^2 299 / 300 ENTER 301 RF 302 RCL 06 303 2 304 + 305 RCL 12 306 * 307 RCL 09 308 RCL 07 309 / 310 SQRT 311 RCL 06 312 * 313 + 314 * 315 2 316 / 317 ST+ 18 318 RCL 15 319 RCL 14 320 RCL 13 321 RF 322 RCL 06 323 * 324 CHS 325 RCL 18 326 + 327 ST+ X 328 STO 04 329 CLX 330 STO 05 331 STO 06         332 STO 07 333 GTO 07 334 LBL 05 335 STO 11 336 RCL 06 337 RCL 08 338 - 339 RCL 06 340 RCL 09 341 - 342 * 343 STO 15 344 RCL 07 345 RCL 11 346 - 347 RCL 07 348 RCL 06 349 - 350 / 351 STO 12 352 RCL 11 353 RCL 08 354 - 355 RCL 06 356 LASTX 357 - 358 FS? 10 359 GTO 00 360 / 361 STO 13 362 RCL 11 363 RCL 09 364 - 365 RCL 06 366 LASTX 367 - 368 / 369 STO 14 370 GTO 02 371 LBL 00 372 * 373 RCL 09 374 X^2 375 + 376 STO 13 377 RCL 11 378 RCL 06 379 ST- Y 380 RCL 08 381 - 382 X^2 383 RCL 09 384 ST* Z 385 X^2 386 + 387 STO 15 388 ST/ 13 389 / 390 STO 14 391 LBL 02 392 1 393 X<>Y 394 RCL 13 395 RCL 12 396 FS? 10 397 RJZ 398 FC? 10 399 RJ 400 RCL 11 401 RCL 06 402 - 403 * 404 3 405 / 406 X<> 14 407 RCL 13 408 RCL 12 409 FS? 10 410 RFZ 411 FC? 10 412 RF 413 RCL 06 414 * 415 RCL 14         416 + 417 RCL 15 418 SQRT 419 / 420 RCL 11 421 RCL 06 422 - 423 RCL 07 424 LASTX 425 - 426 / 427 ABS 428 SQRT 429 ST+ X 430 * 431 RTN 432 LBL 03 433 FC? 09 434 GTO 03 435 RCL 08 436 STO 06 437 RCL 09 438 STO 07 439 RCL 11 440 STO 08 441 RCL 12 442 CHS 443 SQRT 444 STO 09 445 LBL 03 446 FC? 10 447 GTO 03 448 RCL 13 449 STO 08 450 RCL 14 451 CHS 452 SQRT 453 STO 09 454 LBL 03 455 FC? 10 456 FS? 09 457 GTO 04 458 RCL 06 459 RCL 07 460 X>Y? 461 X<>Y 462 RCL 08 463 X>Y? 464 X<>Y 465 RCL 09 466 X>Y? 467 X<>Y 468 STO 06 469 RDN 470 XY 472 RDN 473 X>Y? 474 X<>Y 475 STO 07 476 X<> T 477 XY 479 STO 09 480 X<>Y 481 STO 08 482 1 483 X>Y? 484 GTO 02 485 RCL 07 486 X>Y? 487 GTO 03 488 X<> 06 489 X<> 08 490 STO 07 491 GTO 03 492 LBL 02 493 RCL 08 494 X<> 06 495 STO 08 496 RCL 09 497 X<> 07 498 STO 09 499 GTO 03 500 LBL 04 501 RCL 06 502 RCL 07         503 X>Y? 504 X<>Y 505 STO 06 506 X<>Y 507 STO 07 508 LBL 03 509 FS?C 09 510 SF 10 511 RCL 06 512 1 513 XY? 517 GTO 06 518 X<> 06 519 STO 07 520 1 521 XEQ 05 522 STO 04 523 RCL 06 524 CHS 525 X<=0? 526 GTO 00 527 CLX 528 XEQ 05 529 ST- 04 530 LBL 00 531 CLX 532 STO 05 533 STO 07 534 GTO 07 535 LBL 02 536 CLX 537 XEQ 05 538 STO 04 539 ST+ X 540 STO 05 541 1 542 XEQ 05 543 ST- 04 544 CLX 545 X<> 06 546 STO 07 547 GTO 07 548 LBL 06 549 1 550 XEQ 05 551 STO 04 552 RCL 07 553 XEQ 05 554 ST+ X 555 STO 05 556 RCL 06 557 CHS 558 X<=0? 559 GTO 07 560 CLX 561 XEQ 05 562 ST- 04 563 ST+ X 564 ST- 05 565 LBL 07 566 RCL 06 567 X<0? 568 CLX 569 STO 06 570 DEG 571 90 572 TAN 573 RCL 07 574 X<=0? 575 X<>Y 576 STO 07 577 SF 24 578 RCL 02 579 ABS 580 SQRT 581 ST/ 04 582 ST/ 05 583 RCL 20 584 CHS 585 X#0? 586 SIGN 587 STO 09 588 977.7922214 589 RCL 00         590 / 591 ST* 04 592 ST* 05 593 RCL 20 594 ABS 595 SQRT 596 X=0? 597 SIGN 598 / 599 ST* 06 600 ST* 07 601 STO 08 602 RCL 05 603 RCL 04 604 CF 24 605 END

( 758 bytes / SIZE 021 )

 STACK INPUTS OUTPUTS T H0 > 0 k Z (Omega) mat R0 Y (Omega)lambda P X (Omega)rad T

T = Age of the Universe in Gigayears  ,       k = +1 , 0 , -1  =  spherical, euclidean , hyperbolic Universes
P = Period of the Universe in Gigayears   ,   R0 = current scale factor of the Universe in Gigalightyears

Example1:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001

71    ENTER^
-0.1   ENTER^
-0.003  ENTER^
-0.001  XEQ "AUM"   >>>>   T  = 14.61256009                           ---Execution time = 53s---
RDN   P =  796.1455560
RDN   R0 = 13.10701187
RDN    k = -1 ( hyperbolic space )

-We also have in R06 & R07

Rmin = 1.306380222
Rmax = 250.8400402

Example2:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049   ( our Universe ? )

71      ENTER^
0.044    ENTER^
0.521    ENTER^
0.000049  XEQ "AUM"   >>>>   T  = 15.60221457                            ---Execution time = 54s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 20.88180628
RDN    k = -1

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example3:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,   (Omega) rad  = 0.000049   ( our Universe ? )

71       ENTER^
0.271     ENTER^
0.732     ENTER^
0.000049   XEQ "AUM"   >>>>   T  = 13.66788679                            ---Execution time = 106s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 249.4075046
RDN    k = +1   ( spherical space )

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example4:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.6 ,  (Omega)lambda  = 0.4 ,   (Omega) rad  = 0.1

71     ENTER^
-0.6    ENTER^
0.4    ENTER^
0.1    XEQ "AUM"   >>>>   T  = 17.27498133                             ---Execution time = 44s---
RDN   P =  0   ( no periodic universe )
RDN   R0 = 13.13082118
RDN    k = -1   ( hyperbolic space )

Rmin = 0
Rmax = 9.999999999 E99 = +infinity

Example5:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = 0.01 ,  (Omega)lambda  = -0.4 ,   (Omega) rad  = 0.01

71    ENTER^
0.01   ENTER^
-0.4   ENTER^
0.01   XEQ "AUM"   >>>>   T  = 11.30457458                                ---Execution time = 107s---
RDN   P =  66.21313789
RDN   R0 = 11.72326781
RDN    k = -1 ( hyperbolic space )

Rmin = 0
Rmax = 21.83998798

Example6:      Hubble "constant" = 71 km/s/Mpc  ,  (Omega)mat = -0.4 ,  (Omega)lambda  = +0.001 ,   (Omega) rad  = 1.4

71    ENTER^
-0.4   ENTER^
0.001  ENTER^
1.4    XEQ "AUM"   >>>>   T  = 6.484938348                           ---Execution time = 109s---
RDN   P =  900.0920468
RDN   R0 = 435.5000702
RDN    k = +1 ( spherical space )

Rmin = 0
Rmax = 1889.773264

-In example 6, T is calculated by a difference between 2 numbers close to P
-So, the precision may be only 6 decimals with an HP41:
-With Free42 , it yields:  T = 6.4849377713

6°)
Redshift -> Distance

a)  Gauss-Legendre 2-Point Formula

Data Registers:     R00 = H0                                                 ( Registers R00-R01-R02-R03 are to be initialized before executing "Z-D" )

R01 = (Omega)mat            R04 = D          R07 = R0       R09 to R14:  temp

R02 = (Omega)lambda       R05 = D0         R08 = k
R03 = (Omega)rad            R06 = DL

Flags:  /
Subroutines:  /

 01 LBL "Z-D"  02 DEG  03 STO 14  04 SIGN  05 +  06 STO 10            07 RCL 01  08 LASTX  09 -  10 RCL 02  11 +  12 RCL 03  13 +  14 STO 11  15 X#0?  16 SIGN 17 STO 08  18 977.7922214  19 RCL 00  20 /  21 STO 09  22 RCL 11            23 ABS  24 SQRT  25 X=0?  26 SIGN  27 /  28 STO 07  29 1  30 RCL 10  31 1/X  32 STO 12 33 -  34 RCL 14  35 ST+ 14  36 /  37 STO Y  38 3  39 SQRT  40 /  41 STO 13            42 -  43 STO 06  44 2  45 /  46 ST- 12  47 CLX  48 STO 04 49 STO 05  50 LBL 01  51 RCL 12  52 RCL 13  53 X<> 06  54 STO 13  55 +  56 STO 12            57 ENTER  58 ENTER  59 X^2  60 RCL 02  61 *  62 RCL 11  63 -  64 * 65 RCL 01  66 +  67 *  68 RCL 03  69 +  70 SQRT  71 1/X  72 ST+ 05  73 *  74 ST+ 04  75 DSE 14  76 GTO 01  77 RCL 09            78 RCL 13  79 *  80 60 81 SIN  82 *  83 ST* 04  84 ST* 05  85 RCL 05  86 RCL 07            87 /  88 ENTER  89 E^X-1  90 LASTX  91 CHS  92 E^X-1  93 -  94 2  95 /  96 RCL Y 97 R-D  98 SIN  99 RCL 08           100 X#0? 101 ENTER 102 X>0? 103 ENTER 104 R^ 105 RCL 10 106 * 107 RCL 07 108 * 109 STO 06 110 RCL 05 111 RCL 04 112 END

( 142 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS Z / DL Y z D0 X N D

Example1:      H0 = 68 km/s/Mpc  ,  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049        z = 7

68       STO 00
0.044     STO 01
0.521     STO 02
0.000049  STO 03          and if  you choose   N = 20

7   ENTER^
20  XEQ "Z-D"     >>>>     D = 14.75577029  Gly                             ---Execution time = 42s---
RDN     D0 = 35.73893637  Gly
RDN     DL = 432.2896646  Gly

-With N = 40, it yields:   14.75576636     35.73909040     432.2929577
-With N = 100 -------    14.75576611     35.73910080     432.2931801

-And we have:  R0 = R07 = 21.80306243 Gly  and  k = R08 = -1 ( hyperbolic space )

Example2:      H0 = 68 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001         z = 7

68        STO 00
-0.1       STO 01
-0.003     STO 02
-0.001     STO 03          and if  we choose   N = 20

7   ENTER^
20  XEQ "Z-D"     >>>>     D = 13.88126253  Gly                             ---Execution time = 42s---
RDN     D0 = 35.94478929  Gly
RDN     DL = 752.8766458  Gly

-With N = 40, it yields:   13.88233261     35.95641236     753.5230703
-With N = 100 -------    13.88244594     35.95762620     753.5906104

-And we have:  R0 = R07 = 13.68526238 Gly  and  k = R08 = -1 ( hyperbolic space )

Note:

-The precision is relatively good, provided z is not too close to the maximum redshift.

-The precision is controlled by the display format.

Data Registers:     R00 = H0                                                 ( Registers R00-R01-R02-R03 are to be initialized before executing "Z-D" )

R01 = (Omega)mat            R04 = D          R07 = R0       R09 to R14:  temp

R02 = (Omega)lambda       R05 = D0         R08 = k
R03 = (Omega)rad            R06 = DL

Flags:  /
Subroutines:  /

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

( 235 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS Z / DL Y / D0 X z D

Example1:      H0 = 68 km/s/Mpc  ,  (Omega)mat = 0.044 ,  (Omega)lambda  = 0.521 ,   (Omega) rad  = 0.000049        z = 7

68       STO 00
0.044     STO 01
0.521     STO 02
0.000049  STO 03      and with  FIX 7

7   XEQ "Z-D"     >>>>     D = 14.75576607  Gly                             ---Execution time = 142s---
RDN     D0 = 35.73910110  Gly
RDN     DL = 432.2931865  Gly

-And we have:  R0 = R07 = 21.80306243 Gly  and  k = R08 = -1 ( hyperbolic space )

Example2:      H0 = 68 km/s/Mpc  ,  (Omega)mat = -0.1 ,  (Omega)lambda  = -0.003 ,   (Omega) rad  = -0.001         z = 7

68        STO 00
-0.1       STO 01
-0.003     STO 02
-0.001     STO 03       FIX 7

7   XEQ "Z-D"     >>>>     D = 13.88244965  Gly                             ---Execution time = 145s---
RDN     D0 = 35.95766592  Gly
RDN     DL = 753.5928202  Gly

-And   R0 = R07 = 13.68526238 Gly  and  k = R08 = -1 ( hyperbolic space )

c)  Cosmological Constant = 0 , Positive Matter, Negative Pressure

-These universes have spherical spaces ( k = +1 )
-All the distances are computed with elementary functions.

-You choose the minimum, current & maximum scale factors in R01 to R03

Data Registers:       R00-R07-R08: temp                                               ( Registers R01-R02-R03 are to be initialized before executing "Z-D" )

R01 = Rmin            R04 = D

R02 = R0               R05 = D0
R03 = Rmax            R06 = DL
Flags:  /
Subroutines:  /

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

( 76 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS Z / DL Y / D0 X z D

Example:         Rmin = 2   R0 = 41   Rmax = 257  ( Gygalightyears )       z = 7

2  STO 01
41   STO 02
257  STO 03

7   XEQ "Z-D"   >>>>     D = 11.60902212   Gly                         ---Execution time = 3s---
RDN     D0 = 23.85152165  Gly
RDN     DL = 180.2301790  Gly

d)  No Matter, Negative Pressure

-These universes have hyperbolic spaces ( k = -1 )
-The light-travel time distance is computed with elementary functions.
-The comoving distance employs Carlson elliptic integrals RF

-You choose the minimum, current & maximum scale factors in R01 to R03

Data Registers:       R00 thru R11: temp                                               ( Registers R01-R02-R03 are to be initialized before executing "Z-D" )

R01 = Rmin            R04 = D

R02 = R0               R05 = D0
R03 = Rmax            R06 = DL

Flags:  /
Subroutine:  ( M-Code )  RF  ( cf "Carlson Elliptic Integrals for the HP41" )

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

( 162 bytes / SIZE 012 )

 STACK INPUTS OUTPUTS Z / DL Y / D0 X z D

Example:         Rmin = 0.7   R0 = 14   Rmax = 314  ( Gygalightyears )       z = 7    ,    z = 19 = zmax

0.7  STO 01
14   STO 02
314  STO 03

7   XEQ "Z-D"   >>>>     D = 12.38326745   Gly                                            ---Execution time = 11s---
RDN     D0 = 29.70733010  Gly
RDN     DL = 460.7466900  Gly

19   XEQ "Z-D"   >>>>     D = 13.98718383   Gly
RDN     D0 = 51.64270084  Gly
RDN     DL = 5595.856028  Gly

-With z = 19,  D has the same value as the age of this universe ( more exactly the time since the last minimum scale factor...

e)  No Pressure, Negative Matter

-These universes have hyperbolic spaces ( k = -1 )
-The light-travel time distance and the comoving distance employ Carlson elliptic integrals RF & RJ

-You choose the minimum, current & maximum scale factors in R01 to R03

Data Registers:       R00 thru R12: temp                                               ( Registers R01-R02-R03 are to be initialized before executing "Z-D" )

R01 = Rmin            R04 = D

R02 = R0               R05 = D0
R03 = Rmax            R06 = DL

Flags:  /
Subroutines:  ( M-Code )  RF & RJ  ( cf "Carlson Elliptic Integrals for the HP41" )

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

( 168 bytes / SIZE 0123)

 STACK INPUTS OUTPUTS Z / DL Y / D0 X z D

Example:         Rmin = 1   R0 = 14   Rmax = 314  ( Gygalightyears )   ,    z = 7

1    STO 01
14   STO 02
314  STO 03

7   XEQ "Z-D"   >>>>     D =  13.56149794   Gly                                                ---Execution time = 21s---
RDN     D0 = 33.91425300  Gly
RDN     DL = 626.3433837  Gly

References:

[1]  B. C. Carlson - "A Table of Elliptic Integrals of the Third Kind"
[2]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[3]  B. C. Carlson - "A Table of Elliptic Integrals: two quadratic factors"