hp41programs

Template

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


Overview
 

1°) Empty Universes

  a)  0 < lambda < 1
  b)  General Case

2°) Tolman Universes

  a)  General Case
  b)  Cyclic Universes

3°) Our Universe ?

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

4°)  More General Cyclic Universes  ( 3 programs )

5°)  More General Universes

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

6°)  Redshifht -> Distance

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



Latest Update:   paragraph 6°)


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

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

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

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

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


-In paragraph 6, we calculate:

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

Formulae:

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

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

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


1°) Empty Universes  


 
    a)  0 < Lambda < 1
 

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


Data Registers: /

Flags: /
Subroutines: /


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

 
( 49 bytes / SIZE 000 )

 

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


Example:
     (Omega)lambda = 0.41

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

Note:

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


      b)  General Case



Data Registers: /

Flags: /
Subroutines: /


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

     
         ( 151 bytes / SIZE 000 )

 

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

    Where  P  is the period

Example1:
    

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

Example2:     

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


Example3:     

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


Example4:     

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


Example5:     

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

2°) Tolman Universes    


      a)  General Case


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


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

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


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


         ( 279 bytes / SIZE 011 )

              

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


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

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

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

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

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

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


Notes:

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


      b)  Cyclic Universes


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


Data Registers:       R00 = L (Gy)-2                       

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



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

     
         ( 88 bytes / SIZE 008 )

 

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

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

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

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

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

Note:

 k is always equal to -1  ( hyperbolic space )

3°) Our Universe ?
 

    a) Via Carlson Integrals


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

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

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

-These conditions are precisely satisfied by our Universe !




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


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

       
        ( 230 bytes / SIZE 022 )

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


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

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

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

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

Note:

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


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


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

       
        ( 203 bytes / SIZE 016 )

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


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

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

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

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


   b) Via Numerical Integration  ( 3 programs )
 

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


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

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


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


      ( 118 bytes / SIZE 011 )

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

    Where N = number of subintervals 

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

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

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

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

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

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

Notes:

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

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


Data Registers:     R00 = H                               R04 thru R10: temp         

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


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


      ( 119 bytes / SIZE 011 )

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


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

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

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

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


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

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

Notes:

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

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


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


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

     
         ( 107 bytes / SIZE 008 )

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


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

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

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

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

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

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

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

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

Notes:

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

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

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


4°)  More General Cyclic Universes


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

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

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


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

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


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


         ( 370 bytes / SIZE 019 )

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

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

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

            71      STO 00

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

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

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


Notes:

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

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


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

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


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


         ( 288 bytes / SIZE 016 )

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

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

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

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

  and  R00 = L = -0.000030330 Gy-2


-We can use M-Code routines  RF  &  RJ

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

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

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

  and with this formula, the program becomes slightly shorter.


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

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


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


      ( 148 bytes / SIZE 012 )

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

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

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

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

  and  R00 = L = -0.000030330 Gy-2

Note:

-Cf also "Cyclic Universes for the HP41"


4°) More General Universes

      a)  Cosmological Constant = 0


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



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

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



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

     
         ( 97 bytes / SIZE 007 )

 

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

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

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

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

Note:

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

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



Data Registers:    R00 = H0    

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


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


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


         ( 313 bytes / SIZE 012 )

              

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

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

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

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

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


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

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

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


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

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

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


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

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

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


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

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

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


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

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

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


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

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

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


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

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

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


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

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

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


Notes:

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


      b)  Cosmological Constant # 0



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



Data Registers:    R00 = H0    

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

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

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


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


         ( 649 bytes / SIZE 020+??? )

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

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

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

                      FIX 7

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

-We also have in R06 & R07

    Rmin = 1.306380222
    Rmax = 250.8400402

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

                        FIX 7

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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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

                         FIX 7

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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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

                       FIX 7

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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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

                    FIX 7

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

    Rmin = 0
    Rmax = 21.83998798

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

                    FIX 7

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

    Rmin = 0
    Rmax = 1889.773264




Notes:

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

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

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

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

   30   ENTER^
  100  XEQ 10


Data Registers:    R00 = H0    

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

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

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


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


         ( 631 bytes / SIZE 027 )

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

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

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


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

-We also have in R23 & R24

    Rmin = 1.306380222
    Rmax = 250.8400402

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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

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

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

 ... and so on ...

Note:

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


Data Registers:    R00 = H0    

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

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

-Lines 156-334 are three-byte GTOs


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


         ( 977 bytes / SIZE 021 )

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

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

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


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

-We also have in R06 & R07

    Rmin = 1.306380222
    Rmax = 250.8400402

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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


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

    Rmin = 0
    Rmax = 21.83998798

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


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

    Rmin = 0
    Rmax = 1889.773264


Notes:

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

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

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

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



Data Registers:    R00 = H0    

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

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

-Lines 155-333 are three-byte GTOs


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


         ( 758 bytes / SIZE 021 )

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

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

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


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

-We also have in R06 & R07

    Rmin = 1.306380222
    Rmax = 250.8400402

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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


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

    Rmin = 0
    Rmax = 9.999999999 E99 = +infinity

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


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

    Rmin = 0
    Rmax = 21.83998798

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


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

    Rmin = 0
    Rmax = 1889.773264

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


6°)
Redshift -> Distance  


 
    a)  Gauss-Legendre 2-Point Formula




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

                               
 R01 = (Omega)mat            R04 = D          R07 = R0       R09 to R14:  temp
                               
 R02 = (Omega)lambda       R05 = D0         R08 = k    
                                 R03 = (Omega)rad            R06 = DL

Flags:  /
Subroutines:  /


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


         ( 142 bytes / SIZE 015 )

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

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

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

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

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

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


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

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

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

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

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

Note:

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


      b)  Tanh-Sinh Quadrature  


-The precision is controlled by the display format.


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

                               
 R01 = (Omega)mat            R04 = D          R07 = R0       R09 to R14:  temp
                               
 R02 = (Omega)lambda       R05 = D0         R08 = k    
                                 R03 = (Omega)rad            R06 = DL

Flags:  /
Subroutines:  /


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


         ( 235 bytes / SIZE 017 )

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

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

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

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


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


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

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

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


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


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


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

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


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

                               
 R01 = Rmin            R04 = D       
                               
 R02 = R0               R05 = D0          
                                 R03 = Rmax            R06 = DL
Flags:  /
Subroutines:  /



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

     
         ( 76 bytes / SIZE 009 )

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

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

    2  STO 01
   41   STO 02
  257  STO 03

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



      d)  No Matter, Negative Pressure


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

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


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

                               
 R01 = Rmin            R04 = D       
                               
 R02 = R0               R05 = D0          
                                 R03 = Rmax            R06 = DL

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


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


         ( 162 bytes / SIZE 012 )

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

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

   0.7  STO 01
   14   STO 02
  314  STO 03

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

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

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


      e)  No Pressure, Negative Matter


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

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


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

                               
 R01 = Rmin            R04 = D       
                               
 R02 = R0               R05 = D0          
                                 R03 = Rmax            R06 = DL

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


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


         ( 168 bytes / SIZE 0123)

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

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

    1    STO 01
   14   STO 02
  314  STO 03

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



References:


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