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')² + k.c² = ( -(8.PI.G/c² ).p + L.c² ).R²(t) ( ' = d/dt )
(R')² + k.c² = ( (8.Pi.G/3) (rho) + (L/3).c² ).R²(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 = 10⁹ )

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

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

1°) Radiation Only - Negative Pressure

a) Time <-> Scale Factor

-Here, the density of matter = 0.
-In this case, rho = 3.p_rad/c²

-Einstein's equations become:

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

-Adding these 2 equations yields:

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

-Multiplying by R.R' gives:

2.R².R'.R'' + 2.R.(R')³ - 2.R.R' = 4.(L/3).R³.R' whence

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

( R.R' )² = R² + (L/3).R⁴ + C where C = Constant

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

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

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

if we choose t = 0 when R = R_min

-We also have R_max² = -3/L - R_min²

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 = R_min R04 = R_max
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	/	R_max
Z	/	q
Y	/	H
X	/	R or t

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

Example: L = -0.0000314 (Giga-light-year)^-2 , R_min = 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 R_max = 309.0945506 GLY

-Likewise,

R = 14 XEQ "R-T" >>>>    t = 13.93482969
    RDN H = 69.42018199
   RDN q = -0.008045284
   RDN R_max = 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.R_min²

     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 = H₀ ( km/s/Mpc ) ( Registers R00-R01-R02 are to be initialized before executing "ACU" )

• R01 = lam < 0 R03 = R_min R05: temp
• R02 = q R04 = R_max
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: H₀ = 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 = R_min = 2.950955220 Gigalightyears
R04 = R_max = 251.4187655 Gigalightyears

Notes:

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

-In the following variant, you choose R_min , R₀ , R_max and "AUR" returns T , P , H , q & cosmological constant L is stored in R00

Data Registers: R00 = L

                                        R01 = R²_min R04-R05: temp
                                        R02 = R²
    R03 = R²_max
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	R_min	H
Y	R₀	P
X	R_max	T

Example: R_min = 1 R₀ = 14 R_max = 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')² - 1 = L.R² ( if k = -1 )

-Multiplying by R' gives:

2 R.R'.R'' + (R')³ = R' + L.R².R' but 2 R.R'.R'' + (R')³ = (d/dt) (R.(R')²)    whence

R.(R')² = R + (L/3).R³ + C where C is a constant

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

Thus t = §_Rm^R_M Sqrt [ R / { (-L/3) ( R - R_m ).( R_M - R ).( R - R_- ) } ] ( § = integral )

where R_m = minimum scale factor, R_M = maximum scale factor, R_- = the 3rd root of the cubic polynomial ( R + (L/3).R³ + C )

-We have R_M = (-R_m/2) + sqrt [ -3/L - 3 (-R_m/2)² ]
and R_- = (-R_m/2) - sqrt [ -3/L - 3 (-R_m/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 = R_max R05 to R08: temp
• R02 = R_min 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 , R_min = 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 R_min & R_max are known, we can use the following variant:

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

• R01 = R_min R03 to R07: temp
• R02 = R_max
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:    R_min = 1 , R_max = 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')² - 1 = L.R² with a Runge-Kutta method of order 4

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

• R01 = t₀ • R04 = h/2 R06 to R10: temp
• R02 = R₀ • R05 = N
• R03 = R'₀
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 , R_min = 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 = §_Rm^R Sqrt [ R / { (-L/3) ( R - R_m ).( R_M - R ).( R - R_- ) } ] has 2 singularities for R = R_m & R = R_M
-So it's better to make a change of argument:

-We use: R = R_m + ( R_M - R_m ) Sin² u the integral becomes:

t = §₀^u (2/sqrt(-L/3)) Sqrt [ ( R_m + ( R_M - R_m ) Sin² u ) / ( R_m - R_- + ( R_M - R_m ) Sin² u ) ] du

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

• R01 = L R03 = R_max R05 to R13: temp
• R02 = R_min 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 , R_min = 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 = R_max = 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 = R_max 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 §_Rm^R_M Sqrt [ R / { (-L/3) ( R - R_m ).( R_M - 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 = R_max R00 to R09: temp
• R12 = R_min 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 X<Y?
90 X<>Y
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 , R_min = 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 = R_max = 308.3953433 GLY

Notes:

-With R = R13 = R_max 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 R₁ & R₂ and returns the interval of time t(R₂) - t(R₁)
-Gauss-Legendre 3-point formula is used.

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

• R01 = L R03 = R_max R05 to R15: temp R15 = R₁
• R02 = R_min R04 = R_- R10 = t(R₂) - t(R₁)
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	R₁	H
X	R₂	t(R₂) - t(R₁)

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

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

-0.0000314 STO 01
1.4 STO 02

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

RCL 00 10 R/S >>>> t₂₃ = 5.576698372
RDN H = 90.62167575
RDN q = -0.080274055

RCL 00 15 R/S >>>> t₃₄ = 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 R₂ = R_max

     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 = H₀ ( km/s/Mpc ) ( Registers R00-R01-R02 are to be initialized before executing "ACU" )

• R01 = lam < 0 R03 = R_min R06 to R16: temp
• R02 = q R04 = R_max
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: H₀ = 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 = R_min = 1.209358576 Gigalightyears
R04 = R_max = 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 H₀, 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 = R_min (Gly) R04 = H₀    R06 to R14: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

1 STO 01
14 STO 02
314 STO 03

-If we choose N₁ = 20 & N₂ = 60

20 ENTER^
60 XEQ "ACU" >>>> T =15.49057987 Gy ---Execution time = 3m00s---
RDN P = 993.8651290 Gy
RDN H₀ = 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 = R_min (Gly) R04 = H₀    R06 to R12: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period of the Universe )

Example: R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

1 STO 01
14 STO 02
314 STO 03

-If we choose N₁ = 15 & N₂ = 50

20 ENTER^
60 XEQ "ACU" >>>> T =15.49057987 Gy ---Execution time = 2m25s---
RDN P = 993.8651296 Gy
RDN H₀ = 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 = R_min (Gly) R04 = H₀    R06 to R16: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

1 STO 01
14 STO 02
314 STO 03

-If we choose N₁ = 8 & N₂ = 24

8 ENTER^
24 XEQ "ACU" >>>> T =15.49057986 Gy ---Execution time = 1m45s---
RDN P = 993.8651310 Gy
RDN H₀ = 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 = R_min (Gly) R04 = H₀    R06 to R??: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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: R_min = 1 Gly , R₀ = 14 Gly , R_max = 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 H₀ = 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 = R_min (Gly) R04 = H₀    R06 to R19: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

1 STO 01
14 STO 02
314 STO 03

-If we choose N₁ = 2 & N₂ = 3

2 ENTER^
3 XEQ "ACU" >>>> T = 15.49057995 Gy ---Execution time = 41s---
RDN P = 993.8651306 Gy
RDN H₀ = 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 = R_min     R04 = H R10 = P
                                        R02 = R    R05 = q
    R03 = R_max
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	R_min	H
Y	R₀	Period
X	R_max	Age

Example:    R_min = 1 Gly R₀ = 14 Gly , R_max = 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 = §_Rmin^Rmax R.dR / f(R)

where f(R) = Sqrt [ ( -L/3 ) ( R - R_min ) ( R_max - R ) ( R - R₃ ) ( R - R₄ ) ]

with Gauss-Chebyshev formula, assuming that R_min , R_max & R₃ are known

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

• R01 = R_min R04 = R₄ R05 to R09: temp
• R02 = R_max
• R03 = R₃
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:    R_min = 1 , R_max = 314 , R₃ = -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:    R_min = 1 , R_max = 314 , R₃ = +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 R₃ <= R_min such that R₄= -R_min-R_max-R₃ <= R_min too.

     b) Age + Period of the Universe

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

where f(R) = Sqrt [ ( -L/3 ) ( R - R_min ) ( R_max - R ) { R² + ( R_min + R_max ) R + ( R²_max + R_max R_min + R²_min + 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 = R_min R04 = H R06 to R16: temp
• R02 = R R05 = q
• R03 = R_max
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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: L = ( -1/32971 ) Gy^-2 , R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

32971 CHS 1/X STO 00

1 STO 01
14 STO 02
314 STO 03 and if we choose N₁ = 10 & N₂ = 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: R² + ( R_min + R_max ) R + ( R²_max + R_max R_min + R²_min + 3/L ) = 0 are not between R_min & R_max
-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 , R_min , R & R_max , we can also give R_min , R , R_max & R₃ where R₃ = the 3rd root of the quartic ( the other one is -R_min-R_max-R₃ 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 = R_min • R04 = R₃ R06 = H R08 to R17: temp
• R02 = R R05 = R₄ R07 = q
• R03 = R_max
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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example1:    R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = -0.1

   1 STO 01
14 STO 02
314 STO 03
-0.1 STO 04 and if we choose N₁ = 10 & N₂ = 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:    R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 0.1

   1 STO 01
14 STO 02
314 STO 03
0.1 STO 04 and if we choose N₁ = 10 & N₂ = 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 R₃ such that R₃ & R₄ = -R_min-R_max-R₃ are <= R_min ( otherwise, there will be probably a DATA ERROR somewhere... )
-With R₃ = -R_min or -R_max , 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 = R_min • R04 = R₃ R06 = H R08 to R14: temp
• R02 = R R05 = R₄ R07 = q
• R03 = R_max
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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example1:    R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = -0.1

   1 STO 01
14 STO 02
314 STO 03
-0.1 STO 04 and if we choose N₁ = 15 & N₂ = 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:    R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 0.1

   1 STO 01
14 STO 02
314 STO 03
0.1 STO 04 and if we choose N₁ = 15 & N₂ = 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

§_a^r x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ] re-written:

§₀^+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 = R_min     R04 = R₃ R06 = H    R08 to R14: temp
                                        R02 = R    R05 = -R₄ R07 = q
    R03 = R_max
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	R₃	q
Z	R_min	H
Y	R₀	Period
X	R_max	Age

Example:    R_min = 0.9 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 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 = R_m + ( R_M - R_m ) Sin² phi

Data Registers: R00 = L

                                      R01 = R_min     R04 = R₃ R06 = H R08 = T    R10 to R??: temp
                                        R02 = R    R05 = -R₄ R07 = q R09 = P
    R03 = R_max
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	R₃	q
Z	R_min	H
Y	R₀	Period
X	R_max	Age

Example:    R_min = 0.9 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 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 = R_min • R14 = R₃ < R_min R01 = H R03 to R16: temp
• R12 = R R15 = R₄ < R_min R02 = q
• R13 = R_max
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:    R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = -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:    R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 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

§_a^r x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ] = §₀^+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 R_F [ (b-r)/(b-a) , (r+c)/(a+c) , (r+d)/(a+d) ] + [(r-a)/3] R_J [ (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 = R_min     R04 = R₃ R06 = H R11 = T
                                        R02 = R    R05 = -R₄ R07 = q R10 = P
    R03 = R_max
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	R₃	q
Z	R_min	H
Y	R₀	Period
X	R_max	Age

Example:    R_min = 0.9 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 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 R₃ , R_min , R & R_max , 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 y⁴ + ( 1 + ¶ - 2.q - 3.lam ) y² + 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 = H₀ ( km/s/Mpc ) ( Registers R00-R01-R02-R03 are to be initialized before executing "ACU" )

• R01 = lam < 0 R04 = R_min R06 to R16: temp
• R02 = q R05 = R_max
• 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 X<Y?
130 X<>Y
131 RDN
132 X<Y?
133 X<>Y
134 R^
135 X<Y?
136 X<>Y
137 STO 05
138 RDN
139 STO 04
140 X<>Y

141 STO 06
142 1
143 RCL 04
144 X<Y?
145 GTO 14
146 RCL 06
147 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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: H₀ = 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 N₁ = 20 & N₂ = 60

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

and R04 = R_min = 1.306380222 Gly
R05 = R_max = 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 = y_m + ( y_M - y_m ) Sin² 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

§_a^r x dx / sqrt [(x-a) (b-x) (x+c) (x+d) ] = §₀^+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 R_F [ (b-r)/(b-a) , (r+c)/(a+c) , (r+d)/(a+d) ] + [(r-a)/3] R_J [ (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 = R_min     R04 = R₃ R06 = H R09 = T
                                        R02 = R    R05 = -R₄ R07 = q R10 = P
    R03 = R_max
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	R₃	q
Z	R_min	H
Y	R₀	Period
X	R_max	Age

Example:    R_min = 0.9 Gly , R₀ = 14 Gly , R_max = 314 Gly , R₃ = 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
R15 = (Omega)_rad = 0.0004953126994

4°) Another Change of Variable

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

-Another similar change of variable may be R = R_m + ( R_M - R_m ) x² ( 2 - x )²
-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 = R_min (Gly) R04 = H₀    R06 to R14: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

1 STO 01
14 STO 02
314 STO 03

-If we choose N₁ = 20 & N₂ = 60

20 ENTER^
60 XEQ "ACU" >>>> T =15.49057974 Gy    ---Execution time = 2m39s---
   RDN P = 993.8651288 Gy
   RDN H₀ = 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 = R_min (Gly) R04 = H₀    R06 to R16: temp
                                      • R02 = R₀ (Gly) R05 = q
  • R03 = R_max (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	N₁	Period
X	N₂	Age

Where N₁ = Number of subintervals to compute the 1st integral ( age of the Universe )
N₂ = Number of subintervals to compute the 2nd integral ( Period - age of the Universe )

Example: R_min = 1 Gly , R₀ = 14 Gly , R_max = 314 Gly

1 STO 01
14 STO 02
314 STO 03

-If we choose N₁ = 8 & N₂ = 24

8 ENTER^
24 XEQ "ACU" >>>> T =15.49057988 Gy ---Execution time = 1m33s---
RDN P = 993.8651294 Gy
RDN H₀ = 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 )

hp41programs

Cyclic Universes for the HP-41