hp41programs

Template

Age of the Universe 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
  b)  Via Numerical Integration

4°)  More General Cyclic Universes ( 2 programs )
5°)  
More General Universes

  a)  Cosmological Constant = 0
  b)  Cosmological Constant # 0



Latest Updates:   paragraphs 5°) and 2°)b) and 4°)


-The age of the universe may be computed by the following integral:   
        t = (c/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     ,      c = speed of light   &   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 the programs below ( except in paragraphs 2°)b) 4°) & 5°) )  1/H = 13.77172143  E9  years  which corresponds to H = 71 km/s/Mpc
-The following routines also compute R = the radius ( or scale factor ) of the Universe


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

     
         ( 90 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


   b) Via Numerical Integration
 

-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 = -1G

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 !


4°) More General Cyclic Universes  ( 2 programs )


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

   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


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 "AUM" 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 "AUM"
 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 "AUM"   >>>>   T = 12.28422138   gigayears                                 ---Execution time = 4s---
                                 RDN   P = 813.6724974   gigayears
                                 RDN  H0 = 53.38731059  km/s/Mpc
                                 RDN  q  =  0.569266382



      b)  Cosmological Constant # 0



-This program employs Romberg integration.
-The precision is controlled by the display format ( lines 372 to 375 ).



Data Registers:     R00 = H0 ( km/s/Mpc )                                                         ( Register R00 is to be initialized before executing "AUM" )

                                   R01 to R?? temp
Flags:  F09-F10
Subroutines: /

-Lines 155-237-246 are three-byte GTOs


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


         ( 575 bytes / SIZE 020+??? )

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

   Where  T = Age of the Universe in Gigayears  ,       k = -1 , 0 , +1  = spherical, euclidean , hyperbolic Universe
                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      STO 00        FIX 7

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

-We also have in R04 & R05

    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      STO 00        FIX 7

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


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

            71      STO 00        FIX 7

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


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

            71      STO 00        FIX 7

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


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

            71      STO 00        FIX 7

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

-We also have in R05:    Rmax = 21.83998798


Notes:

-We use Bairstow method ( lines 23 to 93 ) to factorize the quartic polynomial before finding its roots.
-Romberg method = LBL 10  ( lines 301 to 378 )
-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 !