Template

# Cyclic Universes for the HP-41

Overview

1°)  Radiation Only - Negative Pressure

a)  Time <-> Scale Factor
b)  Age + Period of the Universe
( 2 programs )

2°)  Zero-Pressure & Negative density of Matter

a)  Period of the Universe
b)  Time -> Scale Factor
c)  Scale Factor -> Time ( with Gauss-Legendre & Carlson Elliptic Integrals )
d) 2 Scale-factors -> Delta-T
e)  Age + Period of the Universe  ( 7 programs )

3°)  More General Universes

a)  Period of the Universe
b)  Age + Period of the Universe
( 8 programs )
c)  My favorite Focal Program

4°)  Another Change of Variable

Latest Update:   paragraph 3°)c)

-These programs solve Einstein's equations

2 R.R'' + (R')2 + k.c2 = ( -(8.PI.G/c2 ).p + L.c2 ).R2(t)              ( ' = d/dt )
(R')2 + k.c2  = ( (8.Pi.G/3) (rho) + (L/3).c2 ).R2(t)

where  R = scale factor of the universe,  L = cosmological constant,  k = curvature,  G = gravitational constant,  c = speed of light,  rho = density.

-In all the following programs, k = -1 ( hyperbolic space )  and  L < 0

-The units are choosen so that c = 1
-Unit of time = 1 Giga-Year,  unit of distance = 1 Giga-Light-Year    ( Giga = 109 )

-We also have:    Hubble "constant"  H = R'/R  &  deceleration parameter  q = - R.R'' / R'2

R

|                       *
|             *                 *                         *
|         *                        *                *
|       *                            *          *
|    *                                  *     *
|*                                         *
|--------------------------------P---------------------------- t
0

1°)  Radiation Only - Negative Pressure

a) Time <-> Scale Factor

-Here, the density of matter = 0.
-In this case, rho = 3.prad/c2

-Einstein's equations become:

2 R.R'' + (R')2 - 1 = - 8(Pi).G.p + L
(R')2 - 1 =  8(Pi).G.p + (L/3).R2

2.R.R'' + 2.(R')2 - 2 = 4.(L/3).R2

-Multiplying by  R.R'  gives:

2.R2.R'.R'' + 2.R.(R')3 - 2.R.R' = 4.(L/3).R3.R'    whence

d/dt ( R.R' )2  =  2.R.R' + 4.(L/3).R3.R'    whence

( R.R' )2  =  R2 + (L/3).R4 + C           where C = Constant

-Thus     R.dR/dt = SQRT [ (L/3).R4 + R2 + C ]

-So, the scale factor R(t)  may be expressed in terms of elementary functions:

R = Sqrt [ -3/(2.L) + ( Rmin2 + 3/(2.L) )  Cos ( 2.t sqrt(-L/3) ) ]                                 and the inverse function too:
t  =  ( 1 / sqrt (-4.L/3) )  Arc Cos [ ( R2 + 3/(2L) ) / ( Rmin2 + 3/(2L) ) ]

if we choose  t = 0  when  R = Rmin

-We also have  Rmax2 = -3/L - Rmin2

Data Registers:              R00 = P                ( Registers R01 & R02 are to be initialized before executing "T-R" & "R-T" )

•  R01 = L                R03 = R            R05-R06: temp
•  R02 = Rmin           R04 = Rmax
Flags: /
Subroutines: /

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

( 127 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS T / Rmax Z / q Y / H X / R or  t

Where   Rmax = maximum scale factor, q = deceleration parameter ,  H = Hubble "constant"
R = scale factor & t = time since the minimum scale factor  Rmin

Example:
L = -0.0000314 (Giga-light-year)-2  ,    Rmin =  1.4  Giga-light-year

-0.0000314   STO 01
1.4          STO 02

t = 16   XEQ 'T-R"   >>>>    R = 16.05368930     GLY
RDN    H =  60.59307218     km/s/Mpc
RDN    q  = -0.004958575    ( unitless )
RDN   Rmax  =  309.0945506    GLY

-Likewise,

R =  14   XEQ "R-T"   >>>>    t   = 13.93482969
RDN    H =  69.42018199
RDN    q  = -0.008045284
RDN  Rmax  =  309.0945506

-We also have the period in  R00 = P = 971.0591303  Giga-Years

Notes:

-Lines 90-91 give the Hubble "constant" in kilometers by second by megaparsec
-Otherwise, H would be expressed in (Giga-years)-1

-We must choose  -3/L > 2.Rmin2

b) Age + Period of the Universe

-This program takes the current Hubble "constant", cosmological parameter lam and deceleration parameter q in R00  R01 & R02
and returns the age of the Universe, the period and the current scale factor in registers X Y Z.

Data Registers:           •  R00 = H0 ( km/s/Mpc )                               ( Registers R00-R01-R02 are to be initialized before executing "ACU" )

•  R01 = lam < 0         R03 = Rmin                R05: temp
•  R02 = q                   R04 = Rmax
Flags: /
Subroutines: /

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

( 96 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Z / Scale Factor Y / Period X / Age

Example:
H0 = 71 km/s/Mpc   lam = -0.003   q = -0.048

71     STO 00
-0.003  STO 01
-0.048  STO 02      XEQ "ACU"   >>>>       Age of the Universe   = 13.09340736  Gigayears
RDN      Period of the Universe = 789.9097505  Gigayears
RDN       Current Scale Factor  =  13.41429720  Gigalightyears

-We also have  R03 = Rmin = 2.950955220  Gigalightyears
R04 = Rmax = 251.4187655  Gigalightyears

Notes:

-Lines 03 to 26 solve a biquadratic equation:   lam  y4  + ( 1 - q - 2 lam ) y2 + ( q + lam ) = 0
-We must have:  1 - q - 2 lam  > 0

-In the following variant, you choose Rmin , R0 , Rmax  and  "AUR"  returns  T , P , H , q  &  cosmological constant L  is stored in R00

Data Registers:             R00 = L

R01 = R2min           R04-R05: temp
R02 = R2
R03 = R2max
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 ST- T  10 X<>Y  11 X^2  12 ST+ 00 13 STO 01  14 -  15 STO 05  16 R^  17 ST* 05  18 -  19 STO 04  20 RCL 01            21 RCL 03  22 -  23 /  24 ACOS 25 D-R  26 RCL 00  27 SQRT  28 *  29 2  30 /  31 X<> 04  32 RCL 05  33 RCL 00            34 /  35 SQRT  36 RCL 02 37 /  38 STO 05  39 X^2  40 RCL 00            41 *  42 RCL 02  43 *  44 ST+ Y  45 /  46 RCL 05  47 977.7922214 48 *  49 3  50 CHS  51 RCL 00            52 /  53 X<> 00  54 SQRT  55 PI  56 *  57 RCL 04  58 END

( 82 bytes / SIZE 008 )

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

Example:
Rmin = 1    R0 = 14   Rmax = 314

1    ENTER^
14   ENTER^
314  XEQ "AUR"  >>>>     T = 13.96898881  Gy                                                                  ---Execution time = 3s---
RDN      P = 986.4650960  Gy
RDN      H = 69.59427440  km/s/Mpc
RDN      q =  -0.003136335

and   R00 =  L =  -0.000030427  (Gy)-2

2°)  Zero-Pressure & Negative Density of Matter

a) Period of the Universe

-If there is no pressure, Einstein's equation becomes:

2 R.R'' + (R')2 - 1 =  L.R2                ( if k = -1 )

-Multiplying by R'  gives:

2 R.R'.R'' + (R')3  =  R'  + L.R2.R'       but  2 R.R'.R'' + (R')3  =  (d/dt) (R.(R')2)    whence

R.(R')2 = R + (L/3).R3 + C         where  C  is a constant

dR/dt = Sqrt [ ( R + (L/3).R3 + C ) / R ]

Thus   t =  §RmRM  Sqrt [ R / { (-L/3) ( R - Rm ).( RM - R ).( R - R- ) } ]            ( § = integral )

where  Rm = minimum scale factor,  RM = maximum scale factor,  R- = the 3rd root of the cubic polynomial   ( R + (L/3).R3 + C )

-We have   RM = (-Rm/2) + sqrt [ -3/L - 3 (-Rm/2)2 ]
and       R-  =  (-Rm/2) - sqrt [ -3/L - 3 (-Rm/2)2 ]  < 0

-Unfortunately, the integral cannot be expressed in terms of elemntary functions: it's an elliptic integral
-But we can use numerical integration formula to compute these results.

-The following routine employs Gauss-Chebyshev formula:

Data Registers:              R00 = P                ( Registers R01 & R02 are to be initialized before executing "P" )

•  R01 = L                R03 = Rmax            R05 to R08: temp
•  R02 = Rmin           R04 = R-
Flags: /
Subroutines: /

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

( 91 bytes / SIZE 008 )

 STACK INPUT OUTPUT X / P

Example:
L = -0.0000314 (GLY)-2  ,    Rmin =  1.4  GLY

-0.0000314   STO 01
1.4          STO 02

XEQ "P"    >>>>    P = 978.7265794  GY                 ---Execution time = 63s---

Notes:

-Gauss-Chebyshev formula is applied with N = 50 subintervals ( line 31 )
-To apply it with, for instance, N = 100:

100  XEQ 10  >>>>    P = 978.7265794  GY

-If  Rmin & Rmax are known, we can use the following variant:

Data Registers:              R00 = L                ( Registers R01 & R02 are to be initialized before executing "P" )

•  R01 = Rmin           R03 to R07: temp
•  R02 = Rmax
Flags: /
Subroutines: /

 01 LBL "P"  02 STO 05        03 STO 06  04 3  05 RCL 01  06 RCL 02  07 +  08 STO 03  09 STO 07  10 LASTX  11 *  12 RCL 01 13 X^2  14 +  15 /  16 CHS  17 STO 00        18 RCL 02  19 RCL 01  20 -  21 2  22 ST/ 07  23 /  24 STO 04 25 CLX  26 RAD  27 LBL 01  28 RCL 05        29 ST+ X  30 PI  31 ST* Y  32 -  33 RCL 06  34 /  35 COS  36 RCL 04 37 *  38 RCL 07  39 +  40 RCL X  41 RCL 03        42 +  43 /  44 SQRT  45 +  46 DSE 05  47 GTO 01  48 PI 49 *  50 RCL 06  51 /  52 12  53 RCL 00        54 CHS  55 /  56 SQRT  57 *  58 DEG  59 END

( 72 bytes / SIZE 008 )

 STACK INPUT OUTPUT X N P

Where  N = number of subintervals for Gauss-Chebyshev formula

Example:
Rmin = 1 , Rmax = 314

1    STO 01
314  STO 02

-If you choose  N = 120

120   XEQ "P"   >>>>   P = 993.8651304                                       ---Execution time = 2m20s---

And  R00 = L = -0.00003033029693

Note:

-Exact value = 993.8651296  ( rounded to 7 decimals )

b) Time -> Scale Factor

-The following routine "T-R" solves the differential equation   2 R.R'' + (R')2 - 1 =  L.R2    with a Runge-Kutta method of order 4

Data Registers:           •  R00 = L               ( Registers R01 & R02 are to be initialized before executing "T-R" )

•  R01 = t0            •  R04 = h/2            R06 to R10: temp
•  R02 = R0          •  R05 = N
•  R03 = R'0
Flags: /
Subroutines: /

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

( 144 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS Z / q Y / H X / R(t)

Example:
L = -0.0000314 (GLY)-2  ,    Rmin =  1.4  GLY  = R(0)   ,    R'(0) = 0

-Calculate  R(16)

-0.0000314   STO 00

0           STO 01
1.4          STO 02
0           STO 03

-If we choose  h/2 = 0.1 ,  N = 80

0.1          STO 04
80           STO 05        XEQ "T-R"    >>>>    t = 16 GY                                      ---Execution time = 5m07s---
RDN   R = 14.15775360  GLY
RDN   H = 65.48419681  km/s/Mpc
RDN    q = -0.052661760

-With  h/2 = 0.05  it yields:  R(16) = 14.15775380

c) Scale Factor -> Time

-The inverse function may be computed faster with Gauss-Legendre formulae
-"R-T" below employs 3-point formula.

-The integrand       t =  §RmR  Sqrt [ R / { (-L/3) ( R - Rm ).( RM - R ).( R - R- ) } ]      has 2 singularities for R = Rm  &  R = RM
-So it's better to make a change of argument:

-We use:  R = Rm + ( RM - Rm ) Sin2 u      the integral becomes:

t = §0u  (2/sqrt(-L/3))  Sqrt [ ( Rm + ( RM - Rm ) Sin2 u ) / ( Rm - R- + ( RM - Rm ) Sin2 u ) ]  du

Data Registers:              R00 = R                ( Registers R01 & R02 are to be initialized before executing "R-T" )

•  R01 = L                R03 = Rmax            R05 to R13: temp
•  R02 = Rmin           R04 = R-                R10 = t(R)
Flags: /
Subroutines: /

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

( 176 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS Z / q Y / H X R t(R)

Example:
L = -0.0000314 (GLY)-2  ,    Rmin =  1.4  GLY

-Calculate  t(15)

-0.0000314  STO 01
1.4        STO 02

15    XEQ "R-T"   >>>>   t = 16.88695895  GY = R10                                      ---Execution time = 106s---
RDN   H = 61.98903072 km/s/Mpc
RDN    q = -0.048999274

-We also have:   R03 = Rmax = 308.3953433 GLY

Notes:

-Gauss-Legendre 3-point formula is applied to N = 32 subintervals ( line 36 ).
-If you want to use another value, for instance N = 100, press  100  XEQ 10  >>>>  16.88695903  in the example above.

-With R = R03 = Rmax  there will be a DATAERROR line 130 because H = 0  ( and q = -infinity )
-But  RCL 10  will give the semi-period  P/2 = 489.3632912  whence  P = 978.7265824  GY  ( P = 978.7265768  with N = 100 )
-You can also press SF 25 before XEQ "R-T" to avoid this DATA ERROR.

-Instead of a Gauss-Legendre formula, we can also use Carlson elliptic Integrals:
-The integral  §RmRM  Sqrt [ R / { (-L/3) ( R - Rm ).( RM - R ).( R - R- ) } ]    belongs to what is called  [ 1  -1  -1  -1 ]  in reference [2]

-The program below employs this formula and calls  "RJ2"  listed in "Elliptic Integrals for the HP41, paragraph 2°)c-2)

Data Registers:              R10 = R                ( Registers R01 & R02 are to be initialized before executing "R-T" )

•  R11 = L                R13 = Rmax            R00 to R09: temp
•  R12 = Rmin           R14 = R-                R01 = t(R)
Flags: /
Subroutines:  "RJ2"  ( cf "Elliptic Integrals for the HP41, paragraph 2°)c-2) )

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

( 166 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS Z / q Y / H X R t(R)

Example:
L = -0.0000314 (GLY)-2  ,    Rmin =  1.4  GLY

-Calculate  t(15)

-0.0000314  STO 11
1.4        STO 12

15    XEQ "R-T"   >>>>   t = 16.88695894  GY = R01                                      ---Execution time = 27s---
RDN   H = 61.98903072 km/s/Mpc
RDN    q = -0.048999274

-We also have:   R13 = Rmax = 308.3953433 GLY

Notes:

-With R = R13 = Rmax  there will be a DATAERROR line 134 because H = 0  ( and q = -infinity )
-But  RCL 01  will give the semi-period  P/2 = 489.3632897  whence  P = 978.7265794  GY
-You can also press SF 25 before XEQ "R-T" to avoid this DATA ERROR.

d) 2 Scale Factors -> Delta-T

-"2R-DT"  takes 2 values of the scale factor R1 & R2  and returns the interval of time  t(R2) - t(R1)
-Gauss-Legendre 3-point formula is used.

Data Registers:              R00 = R2                ( Registers R01 & R02 are to be initialized before executing "2R-DT" )

•  R01 = L                R03 = Rmax            R05 to R15: temp         R15 = R1
•  R02 = Rmin           R04 = R-                R10 =  t(R2) - t(R1)
Flags: /
Subroutines: /

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

( 192 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS Z / q Y R1 H X R2 t(R2) - t(R1)

Example:
L = -0.0000314 (GLY)-2  ,    Rmin =  1.4  GLY

-Calculate t(5)  t(10)  &  t(15)

-0.0000314  STO 01
1.4        STO 02

RCL 02   5   XEQ "2R-DT"   >>>>  t12 = 5.994335857                               ---Execution time = 55s---
RDN   H = 165.9073489
RDN    q =  -0.194146039

RCL 00  10   R/S   >>>>    t23 = 5.576698372
RDN    H =   90.62167575
RDN     q =   -0.080274055

RCL 00  15   R/S   >>>>    t34 =  5.315924725
RDN     H =  61.98903072
RDN      q =  -0.048999274

-So,   t(5)  = 5.994335857
t(10) = 5.994335857 + 5.576698372 = 11.57103423
t(15) = 5.994335857 + 5.576698372 + 5.315924725 = 16.88695896

Notes:

-All the times are expressed in Giga-years,
all the distances in Giga-light-years and the cosmological constant in (Giga-light-years)-2

-The Hubble "constant" in km/s/Mpc

-The interval of integration is divided in 16 subintervals ( line 46 )
-If you want to recalculate with another value - for instance 50 - press  50 XEQ 10

-Press SF 25 before XEQ "2R-DT" to avoid a DATA ERROR if  R2 = Rmax

e) Age + Period of the Universe

-This program takes the current Hubble "constant", cosmological parameter lam and deceleration parameter q in R00  R01 & R02
and returns the age of the Universe, the period and the current scale factor in registers X Y Z.

Data Registers:           •  R00 = H0 ( km/s/Mpc )                               ( Registers R00-R01-R02 are to be initialized before executing "ACU" )

•  R01 = lam < 0          R03 = Rmin                R06 to R16: temp
•  R02 = q                   R04 = Rmax
R05 = R-
Flags: /
Subroutines: /

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

( 224 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS Z / Scale Factor Y / Period X / Age

Example:
H0 = 71 km/s/Mpc   lam = -0.003   q = -0.048

71     STO 00
-0.003  STO 01
-0.048  STO 02      XEQ "ACU"   >>>>       Age of the Universe   = 14.73930565  Gigayears                     ---Execution time = 5m15s---
RDN      Period of the Universe = 796.4604549  Gigayears
RDN       Current Scale Factor  =  13.10107975  Gigalightyears

-We also have  R03 = Rmin = 1.209358576  Gigalightyears
R04 = Rmax = 250.8292223  Gigalightyears

Notes:

-Lines 03 to 46 solve a cubic equation.

-The integrals are computed with Gauss-Legendre 3-point formula and 32 & 64 subintervals ( lines 123 &131 )
-If a lower precision is enough, these lines may be replaced by smaller values:

-With 10 & 20  instead of 32 & 64, it yields:  14.73930567 & 796.4604560  ( almost as accurate ! )

-Instead of the previous inputs, we can also choose the minimum & maximum scale factors and the current "radius" of the Universe.
-The variant below returns the age & the period of the Universe and the current Hubble constant H0, the deceleration parameter q & the cosmological constant L

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R14: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: /
Subroutines: /

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

( 167 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

-If we choose  N1 = 20  &  N2 = 60

20   ENTER^
60   XEQ "ACU"   >>>>   T  =15.49057987  Gy                                                            ---Execution time = 3m00s---
RDN   P  = 993.8651290 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.000030330 (Gy)-2

-The exact values are  T = 15.490579844....  &  P = 993.865129627....  ( obtained with free42 and larger N-values )

Notes:

-This version employs 2-point Gauss-Legendre formula: it almost minimizes the number of bytes.
-We can write a slighly shorter routine if we compute the integral giving the period independantly of the integral giving the age:

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R12: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: F07
Subroutines: /

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

( 153 bytes / SIZE 013 )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period of the Universe )

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

-If we choose  N1 = 15  &  N2 = 50

20   ENTER^
60   XEQ "ACU"   >>>>   T  =15.49057987  Gy                                                            ---Execution time = 2m25s---
RDN   P  = 993.8651296 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.000030330 (Gy)-2

Note:

-A better precision with less function-evaluations could be obtained with Gauss-Legendre 3-point formula:

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R16: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: /
Subroutines: /

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

( 189 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

-If we choose  N1 = 8  &  N2 = 24

8    ENTER^
24   XEQ "ACU"   >>>>   T  =15.49057986  Gy                                                            ---Execution time = 1m45s---
RDN   P  = 993.8651310 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.00003033029693 (Gy)-2

Note:

-The following variant employs Romberg integration

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R??: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: /
Subroutines: /

-Line 15 is a three-byte GTO

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

( 230 bytes / SIZE 017+??? )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y / Period X / Age

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

FIX 7       XEQ "ACU"   >>>>   T  =15.49057996  Gy                                                            ---Execution time = 2m21s---
RDN   P  = 993.8651316 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.000030330 (Gy)-2

Notes:

-The HP41 displays the successive approximations.
-The precision is ccontrolled by the display format.
-FIX 9 could produce infinite loops...

-In the following program, we use Gauss-Legendre 6-point formula

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R19: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: /
Subroutines: /

-Line 11 is a three-byte GTO 03

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

( 298 bytes / SIZE 020 )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

-If we choose  N1 = 2  &  N2 = 3

2   ENTER^
3   XEQ "ACU"   >>>>     T  =  15.49057995  Gy                                                            ---Execution time = 41s---
RDN   P  = 993.8651306 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.000030330 (Gy)-2

-The following "AUM" employs a method similar to the program listed in §3°)c)

Data Registers:              R00 = L

R01 = Rmin           R04 = H       R10 = P
R02 = R                R05 = q
R03 = Rmax
Flag:  F10
Subroutines:  /

-Line 150 is a three-byte GTO 02

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

( 266 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T / q Z Rmin H Y R0 Period X Rmax Age

Example:
Rmin = 1 Gly   R0 = 14 Gly  ,   Rmax = 314 Gly

1     ENTER^
14    ENTER^
314   XEQ "AUM"   >>>>    T = 15.49057986  Gy                                               ---Execution time = 63s---
RDN    P = 993.8651299  Gy
RDN    H = 67.22990236 km/s/Mpc
RDN    q = -0.036404801

and  R00 = L = -0.000030330  (Gy)-2

Notes:

-Line 149  X>Y?  may be replaced by the M-Code  X#Y??  ( cf "A few M-Code routines for the HP41" §8°)c) )

3°)  More General Universes

a) Period of the Universe

-We compute  P = §RminRmax  R.dR / f(R)

where  f(R) = Sqrt [ ( -L/3 ) ( R - Rmin ) ( Rmax - R ) ( R - R3 ) ( R - R4 ) ]

with Gauss-Chebyshev formula, assuming that  Rmin , Rmax & R3  are known

Data Registers:              R00 = L                ( Registers R01-R02-R03 are to be initialized before executing "P" )

•  R01 = Rmin           R04 = R4            R05 to R09: temp
•  R02 = Rmax
•  R03 = R3
Flags: /
Subroutines: /

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

( 89 bytes / SIZE 010 )

 STACK INPUT OUTPUT X N P

Where  N = number of subintervals for Gauss-Chebyshev formula

Example1:
Rmin = 1 , Rmax = 314  ,  R3 = -0.1

1    STO 01
314  STO 02
-0.1   STO 03

-If you choose  N = 120

120   XEQ "P"   >>>>   P = 993.0388652                                       ---Execution time = 2m38s---

And  R00 = L = -0.00003033995617

Example2:      Rmin = 1 , Rmax = 314  ,  R3 = +0.1

1    STO 01
314  STO 02
0.1   STO 03

-If you choose  N = 120

120   XEQ "P"   >>>>   P = 994.7216216

And  R00 = L = -0.00003032063771

Note;

-Choose R3 <= Rmin  such that  R= -Rmin-Rmax-R3 <= Rmin  too.

b) Age + Period of the Universe

-The following program computes  t = §  R.dR / f(R)

where  f(R) = Sqrt [ ( -L/3 ) ( R - Rmin ) ( Rmax - R ) { R2 + ( Rmin + Rmax ) R + ( R2max + Rmax Rmin + R2min + 3/L ) } ]

-It also returns the Hubble "constant" H ( in km/s/Mpc ) & the deceleration parameter q

Data Registers:           •  R00 = L < 0          ( Registers R00-R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin           R04 = H                       R06 to R16: temp
•  R02 = R                R05 = q
•  R03 = Rmax
Flags: /
Subroutines: /

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

( 248 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / q Z / H Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
L = ( -1/32971 ) Gy-2  ,   Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

32971  CHS  1/X   STO 00

1     STO 01
14    STO 02
314   STO 03    and if we choose   N1 = 10  &   N2 = 30

10    ENTER^
30    XEQ "ACU"   >>>>    T = 15.50327701 Gy                                               ---Execution time = 2m24s---
RDN     P = 993.9185580 Gy
RDN     H = 67.21462506 km/s/Mpc
RDN     q  = -0.036631247

Notes:

-This program checks that the roots - if any - of the quadratic equation:   R2 + ( Rmin + Rmax ) R + ( R2max + Rmax Rmin + R2min + 3/L ) = 0  are not between  Rmin & Rmax
-If these conditions are not satisfied, the HP41 displays  NONEXISTENT  ( line 37 )

-Gauss-Legendre 3-point formula is employed to estimate the integrals.

-Instead of  L , Rmin , R & Rmax , we can also give  Rmin , R , Rmax & R3   where   R3 = the 3rd root of the quartic ( the other one is  -Rmin-Rmax-R3 since the sum of the 4 roots = 0 )

Data Registers:              R00 = L                 ( Registers R01-R02-R03-R04 are to be initialized before executing "ACU" )

•  R01 = Rmin         •  R04 = R3       R06 = H                    R08 to R17: temp
•  R02 = R                 R05 = R4       R07 = q
•  R03 = Rmax
Flags: /
Subroutines: /

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

( 241 bytes / SIZE 018 )

 STACK INPUTS OUTPUTS T / q Z / H Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example1:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = -0.1

1     STO 01
14    STO 02
314   STO 03
-0.1   STO 04  and if we choose   N1 = 10  &   N2 = 30

10    ENTER^
30    XEQ "ACU"   >>>>    T = 15.29789632 Gy                                               ---Execution time = 2m28s---
RDN     P = 993.0388654 Gy
RDN     H = 67.47006839 km/s/Mpc
RDN     q  = -0.032865171

and  R00 = L = -0.00003033995617  Gy-2

Example2:      Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

1     STO 01
14    STO 02
314   STO 03
0.1    STO 04  and if we choose   N1 = 10  &   N2 = 30

10    ENTER^
30    XEQ "ACU"   >>>>    T = 15.69777974 Gy                                               ---Execution time = 2m28s---
RDN     P = 994.7216224 Gy
RDN     H = 66.98887584 km/s/Mpc
RDN     q  = -0.039995457

and  R00 = L = -0.00003032063771  Gy-2

Notes:

-Choose  R3 such that  R3 & R4 = -Rmin-Rmax-R3 are <= Rmin   ( otherwise, there will be probably a DATA ERROR somewhere... )
-With  R3 = -Rmin or -Rmax , you will have Universes of paragraph 1

-Register R17 may be replaced by synthetic register M.

-Gauss-Legendre 2-point formula will give a shorter routine:

Data Registers:              R00 = L                 ( Registers R01-R02-R03-R04 are to be initialized before executing "ACU" )

•  R01 = Rmin         •  R04 = R3       R06 = H                    R08 to R14: temp
•  R02 = R                 R05 = R4       R07 = q
•  R03 = Rmax
Flag:  F07
Subroutines: /

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

( 197 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS T / q Z / H Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example1:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = -0.1

1     STO 01
14    STO 02
314   STO 03
-0.1   STO 04  and if we choose   N1 = 15  &   N2 = 60

15    ENTER^
60    XEQ "ACU"   >>>>    T = 15.29789634 Gy                                               ---Execution time = 3m12s---
RDN     P = 993.0388645 Gy
RDN     H = 67.47006839 km/s/Mpc
RDN     q  = -0.032865171

and  R00 = L = -0.00003033995617  Gy-2

Example2:      Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

1     STO 01
14    STO 02
314   STO 03
0.1    STO 04  and if we choose   N1 = 15  &   N2 = 60

15    ENTER^
60    XEQ "ACU"   >>>>    T = 15.69777977 Gy
RDN     P = 994.7216227 Gy
RDN     H = 66.98887585 km/s/Mpc
RDN     q  = -0.039995458

and  R00 = L = -0.00003032063771  Gy-2

Note:

-We can also use Exp-Sinh quadrature to compute the integral

§ar x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ]       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)

Data Registers:              R00 = L

R01 = Rmin           R04 = R3         R06 = H                    R08 to R14: temp
R02 = R                R05 = -R4       R07 = q
R03 = Rmax
Flags: /
Subroutines: /

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

( 270 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T R3 q Z Rmin H Y R0 Period X Rmax Age

Example:
Rmin = 0.9 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

SCI 6

0.1   ENTER^
0.9   ENTER^
14    ENTER^
314   XEQ "ACU"   >>>>    T = 15.60553775  Gy                                               ---Execution time = 7mn00s---
RDN    P = 994.0878162  Gy
RDN    H = 67.24655367 km/s/Mpc
RDN    q =  -0.035891530

and  R00 = L = -0.000030330  (Gy)-2

Notes:

-The precision is controlled by the display setting.
-Neither the shortest nor the fastest program !

-Here is a program using Romberg integration and the change of variable:  R = Rm + ( RM - Rm ) Sin2 phi

Data Registers:              R00 = L

R01 = Rmin           R04 = R3         R06 = H     R08 = T               R10 to R??: temp
R02 = R                R05 = -R4       R07 = q      R09 = P
R03 = Rmax
Flags: /
Subroutines: /

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

( 274 bytes / SIZE 018+??? )

 STACK INPUTS OUTPUTS T R3 q Z Rmin H Y R0 Period X Rmax Age

Example:
Rmin = 0.9 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

SCI 7

0.1   ENTER^
0.9   ENTER^
14    ENTER^
314   XEQ "ACU"   >>>>    T = 15.60553776  Gy                                               ---Execution time = 4mn44s---
RDN    P = 994.0878168  Gy
RDN    H = 67.24655367 km/s/Mpc
RDN    q =  -0.035891530

and  R00 = L = -0.000030330  (Gy)-2

Notes:

-The precision is controlled by the display setting.

-The following variant employs Carlson elliptic integrals:

Data Registers:              R00 = L                 ( Registers R11-R12-R13-R14 are to be initialized before executing "ACU" )

•  R11 = Rmin         •  R14 = R3 < Rmin      R01 = H                    R03 to R16: temp
•  R12 = R                 R15 = R4 < Rmin      R02 = q
•  R13 = Rmax
Flags: /
Subroutines:   "RF" & "RJ2"  ( cf "Elliptic Integrals for the HP41" )

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

( 225 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / q Z / H Y / Period X / Age

Example1:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = -0.1

1     STO 11
14    STO 12
314   STO 13
-0.1   STO 14

XEQ "ACU"   >>>>    T = 15.29789630 Gy                                               ---Execution time = 71s---
RDN     P = 993.0388640 Gy
RDN     H = 67.47006839 km/s/Mpc
RDN     q  = -0.032865171

and  R00 = L = -0.00003033995617  Gy-2

Example2:      Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

1     STO 11
14    STO 12
314   STO 13
0.1    STO 14

XEQ "ACU"   >>>>    T = 15.69777975 Gy
RDN     P = 994.7216216 Gy
RDN     H = 66.98887585 km/s/Mpc
RDN     q  = -0.039995458

and  R00 = L = -0.00003032063771  Gy-2

Notes:

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

§ar x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ]  =  §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

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

and with this formula, the program becomes slightly shorter.

Data Registers:              R00 = L

R01 = Rmin           R04 = R3         R06 = H     R11 = T
R02 = R                R05 = -R4       R07 = q      R10 = P
R03 = Rmax
Flag:  F07
Subroutines:  M-Code routines "RF" & "RJ"  ( cf "Carlson Elliptic Integrals for the HP41" )

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

( 194 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T R3 q Z Rmin H Y R0 Period X Rmax Age

Example:
Rmin = 0.9 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

0.1   ENTER^
0.9   ENTER^
14    ENTER^
314   XEQ "ACU"   >>>>    T = 15.60553773  Gy                                               ---Execution time = 42s---
RDN    P = 994.0878143  Gy
RDN    H = 67.24655367 km/s/Mpc
RDN    q = -0.035891530

and  R00 = L = -0.000030330  (Gy)-2

Notes:

-The M-Code routines RF & RJ may be replaced by the focal programs "RF" & "RJ" but in this case, shift the data registers to avoid a register conflict...

-Instead of  L or R3 , Rmin , R & Rmax ,  the following variant takes the Hubble "constant" H, the cosmological parameter lam , the deceleration parameter q & the radiation parameter ¶
-Here, we assume that the quartic equation:   lam y4 + ( 1 + ¶ - 2.q - 3.lam ) y2 + 2 ( q + lam - ¶ ) y + ¶ = 0   has 4 real roots

-The radiation parameter  ¶ # 0  ,  cosmological parameter  lam < 0  and the density of matter   2 ( q + lam - ¶ )  <  0

Data Registers:           •  R00 = H0 ( km/s/Mpc )                               ( Registers R00-R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = lam < 0          R04 = Rmin                 R06 to R16: temp
•  R02 = q                   R05 = Rmax
•  R03 = ¶
Flags: /
Subroutines: /

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

( 302 bytes / SIZE 019 )

 STACK INPUTS OUTPUTS Z / Current Radius Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
H0 = 71 km/s/Mpc  ,  lam = -0.003  ,   q = -0.048  ,   ¶ = -0.001

71   STO 00

-0.003   STO 01
-0.048   STO 02
-0.001   STO 03

-If we choose  N1 = 20  &  N2 = 60

20   ENTER^
60   XEQ "ACU"   >>>>   T  = 14.61256009  Gy                                                            ---Execution time = 3m26s---
RDN   P  = 796.1455560 Gy
RDN  R0 = 13.10701187  Gly

and    R04 = Rmin = 1.306380222  Gly
R05 = Rmax = 250.8400402  Gly

Notes:

-The first 2 roots of the quartic equation are found by Newton's method.
-There is perhaps a risk of infinte loop...
-The integrals arre computed with Gauss-Legendre 2-point formula and  the change of argument  y = ym + ( yM - ym ) Sin2 u

-The results obtained in the example above are very close to those obtained in paragraph 2°e) where ¶ = 0

-Cf also "Age of the Universe for the HP41" paragraph 4°)

c) My Favorite Focal Program

-This program also employs the change of variable:    x = (a.y+r) / (y+1)   which yields

§ar x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ]  =  §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

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

-"RF" & "RJ" are computed simultaneously: only 1 loop ( lines 47 to 153 ) to calculate  {  ....   }  above

Data Registers:              R00 = L

R01 = Rmin           R04 = R3         R06 = H     R09 = T
R02 = R                R05 = -R4       R07 = q      R10 = P
R03 = Rmax
Flag:  F10
Subroutines:  /

-Line 153 is a three-byte GTO 02

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

( 308 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T R3 q Z Rmin H Y R0 Period X Rmax Age

Example:
Rmin = 0.9 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1

0.1   ENTER^
0.9   ENTER^
14    ENTER^
314   XEQ "AUM"   >>>>    T = 15.60553774  Gy                                               ---Execution time = 63s---
RDN    P = 994.0878149  Gy
RDN    H = 67.24655367 km/s/Mpc
RDN    q = -0.035891530

and  R00 = L = -0.000030330  (Gy)-2

Notes:

-Line 152  X>Y?  may be replaced by the M-Code  X#Y??  ( cf "A few M-Code routines for the HP41" §8°)c) )
-Execution time becomes 61 seconds in example above.

-We can add the following lines after line 222 ( STO 12 ) to compute the current parameters:   (Omega)mat ,  (Omega) lambda & (Omega)rad   in registers R13-R14-R15  respectively:

 223 RCL 01 224 RCL 03           225 * 226 STO 14 227 RCL 04 228 * 229 STO 13 230 RCL 05 231 * 232 CHS 233 RCL 02           234 X^2 235 X^2 236 / 237 STO 15 238 RCL 04 239 RCL 06 240 * 241 RCL 14           242 + 243 RCL 05 244 * 245 RCL 13 246 - 247 RCL 02 248 ST/ Y 249 X^2 250 / 251 RCL 00           252 RCL 12 253 X^2 254 3 255 * 256 / 257 STO 14           258 ST* 15 259 * 260 STO 13

-With this modification, the program has now  290 lines & 348 bytes.
-In the example above, it yields:

R13 = (Omega)mat       =  -0.07704871226
R14 = (Omega) lambda =  -0.002137513639

4°)  Another Change of Variable

-In several programs, I've used the change of variable:  R = Rm + ( RM - Rm ) Sin2 u
-But trigonometric functions are very slow on an HP41.

-Another similar change of variable may be  R = Rm + ( RM - Rm ) x2 ( 2 - x )2
-The results have approximately the same precision and the routines are slightly faster:

-Here are 2 variants of programs listed in paragraph 2°-e)

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R14: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: /
Subroutines: /

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

( 184 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

-If we choose  N1 = 20  &  N2 = 60

20   ENTER^
60   XEQ "ACU"   >>>>   T  =15.49057974  Gy                                                            ---Execution time = 2m39s---
RDN   P  = 993.8651288 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.00003033029693 (Gy)-2

Notes:

-Execution time = 2m39s instead of 3m00s
-This version employs Gauss-Legandre 2-point formula.

-In the following one, Gauss-Legendre 3-point formula is used:

Data Registers:              R00 = L (Gy)-2                                    ( Registers R01-R02-R03 are to be initialized before executing "ACU" )

•  R01 = Rmin (Gly)           R04 = H0       R06 to R16: temp
•  R02 = R0    (Gly)           R05 = q
•  R03 = Rmax (Gly)
Flags: /
Subroutines: /

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

( 206 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / Deceleration parameter Z / Hubble Constant Y N1 Period X N2 Age

Where   N1 = Number of subintervals to compute the 1st integral ( age of the Universe )
N2 = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example:
Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly

1    STO 01
14   STO 02
314  STO 03

-If we choose  N1 = 8  &  N2 = 24

8   ENTER^
24  XEQ "ACU"   >>>>   T  =15.49057988  Gy                                                            ---Execution time = 1m33s---
RDN   P  = 993.8651294 Gy
RDN  H0 =  67.22990237 km/s/Mpc = R04
RDN   q  =  -0.036404801                 = R05

and  R00 =  L = -0.00003033029693 (Gy)-2

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  B. C. Carlson - A Table of Elliptic Integrals of the Third Kind - MATHEMATICS OF COMPUTATION, VOLUME 51, NUMBER 183 JULY 1988, PAGES 267-280
[3]  J. S. Farnes - A unifying theory of dark energy and dark matter: Negative masses and matter creation within a modified CDM framework - Astronomy & Astrophysic  620, A92 (2018)
[4]  Jacques Colin, Roya Mohayaee, Mohamed Rameez, and Subir Sarkar - "Evidence for anisotropy of cosmic acceleration" - Astronomy & Astrophysic 631, L13 (2019)
[5]  Helge Kragh - Cyclic models of the relativistic universe: the early history.
[6]  Jean-E. Charon - Complex Relativity.

( Reference [6] is a controversial unitary theory, but with many interesting ideas )