hp41programs

Template

Cyclic Universes for the HP-41


Overview
 

 1°)  Radiation Only - Negative Pressure

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

 2°)  Zero-Pressure & Negative density of Matter

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

 3°)  More General Universes


   a)  Period of the Universe
   b)  Age + Period of the Universe
 ( 5 programs )

 4°)  Another Change of Variable


-These programs solve Einstein's equations

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

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

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

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

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


 R

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

1°)  Radiation Only - Negative Pressure
 

       a) Time <-> Scale Factor


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

-Einstein's equations become:

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

-Adding these 2 equations yields:

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

-Multiplying by  R.R'  gives:

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

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

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

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

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

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

 if we choose  t = 0  when  R = Rmin

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


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

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

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

 
      ( 127 bytes / SIZE 007 )

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

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

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

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

-Likewise,

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

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

Notes:

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

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


       b) Age + Period of the Universe


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


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

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


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

 
      ( 96 bytes / SIZE 006 )

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

 
Example:
     H0 = 71 km/s/Mpc   lam = -0.003   q = -0.048
 
        71     STO 00
    -0.003  STO 01
    -0.048  STO 02      XEQ "ACU"   >>>>       Age of the Universe   = 13.09340736  Gigayears      
                                                          RDN      Period of the Universe = 789.9097505  Gigayears
                                                          RDN       Current Scale Factor  =  13.41429720  Gigalightyears

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

Notes:

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

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


Data Registers:             R00 = L 

                                        R01 = R2min           R04 = T         R06 = H         
                                        R02 = R2                R05 = P         R07 = q
                                        R03 = R2max
Flags: /
Subroutines: /

 
 01 LBL "ACU"
 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 007 )

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

 
Example:
 

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

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


2°)  Zero-Pressure & Negative Density of Matter


       a) Period of the Universe


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

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

-Multiplying by R'  gives:

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

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

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

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

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

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

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

-The following routine employs Gauss-Chebyshev formula:


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

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


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

 
      ( 91 bytes / SIZE 008 )

 
            STACK             INPUT          OUTPUT
                 X                 /                  P


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

   -0.0000314   STO 01
          1.4          STO 02
           
         XEQ "P"    >>>>    P = 978.7265794  GY                 ---Execution time = 63s---
 
Notes:

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

         100  XEQ 10  >>>>    P = 978.7265794  GY  

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

 

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

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


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

 
      ( 72 bytes / SIZE 008 )

 
            STACK             INPUT          OUTPUT
                 X                 N                  P

     Where  N = number of subintervals for Gauss-Chebyshev formula

Example:
     Rmin = 1 , Rmax = 314
 
    1    STO 01
  314  STO 02

-If you choose  N = 120

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

 And  R00 = L = -0.00003033029693

Note:

-Exact value = 993.8651296  ( rounded to 7 decimals )


       b) Time -> Scale Factor


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




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

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


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

 
      ( 144 bytes / SIZE 011 )

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

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

-Calculate  R(16)

   -0.0000314   STO 00

            0           STO 01
          1.4          STO 02
            0           STO 03

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

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

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


       c) Scale Factor -> Time


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

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

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

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


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

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


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

 
      ( 176 bytes / SIZE 014 )

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

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

-Calculate  t(15)

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

-We also have:   R03 = Rmax = 308.3953433 GLY
   
Notes:  

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

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


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

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


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

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



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

 
      ( 166 bytes / SIZE 015 )

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

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

-Calculate  t(15)

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

-We also have:   R13 = Rmax = 308.3953433 GLY
   
Notes:  

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


       d) 2 Scale Factors -> Delta-T


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



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

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


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

 
      ( 192 bytes / SIZE 016 )

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

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

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

   -0.0000314  STO 01
          1.4        STO 02

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

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

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


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

Notes:

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

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

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

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


       e) Age + Period of the Universe


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


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

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


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

 
      ( 224 bytes / SIZE 017 )

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

 
Example:
     H0 = 71 km/s/Mpc   lam = -0.003   q = -0.048
 
        71     STO 00
    -0.003  STO 01
    -0.048  STO 02      XEQ "ACU"   >>>>       Age of the Universe   = 14.73930565  Gigayears                     ---Execution time = 5m15s---
                                                          RDN      Period of the Universe = 796.4604549  Gigayears
                                                          RDN       Current Scale Factor  =  13.10107975  Gigalightyears

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

Notes:

-Lines 03 to 46 solve a cubic equation.

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

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

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


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

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


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

   
      ( 167 bytes / SIZE 015 )

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

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

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

    1    STO 01
   14   STO 02
  314  STO 03

-If we choose  N1 = 20  &  N2 = 60

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

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

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


Notes:

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


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

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


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

   
      ( 153 bytes / SIZE 013 )

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

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

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

    1    STO 01
   14   STO 02
  314  STO 03

-If we choose  N1 = 15  &  N2 = 50

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

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


Note:

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


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

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


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

   
      ( 189 bytes / SIZE 016 )

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

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

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

    1    STO 01
   14   STO 02
  314  STO 03

-If we choose  N1 = 8  &  N2 = 24

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

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


Note:

-The following variant employs Romberg integration

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

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

-Line 15 is a three-byte GTO

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

 
      ( 230 bytes / SIZE 017+??? )

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

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

    1    STO 01
   14   STO 02
  314  STO 03  

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

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

Notes:

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

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


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

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


-Line 11 is a three-byte GTO 03

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

   
      ( 298 bytes / SIZE 020 )

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

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

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

    1    STO 01
   14   STO 02
  314  STO 03

-If we choose  N1 = 2  &  N2 = 3

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

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


3°)  More General Universes


       a) Period of the Universe


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

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

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


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

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


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

 
      ( 89 bytes / SIZE 010 )

 
            STACK             INPUT          OUTPUT
                 X                 N                  P

     Where  N = number of subintervals for Gauss-Chebyshev formula

Example1:
     Rmin = 1 , Rmax = 314  ,  R3 = -0.1
 
    1    STO 01
  314  STO 02
 -0.1   STO 03

-If you choose  N = 120

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

 And  R00 = L = -0.00003033995617


Example2:      Rmin = 1 , Rmax = 314  ,  R3 = +0.1
 
    1    STO 01
  314  STO 02
  0.1   STO 03

-If you choose  N = 120

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

 And  R00 = L = -0.00003032063771

Note;

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


       b) Age + Period of the Universe


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

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

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


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

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


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

   
      ( 248 bytes / SIZE 017 )

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

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

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

    32971  CHS  1/X   STO 00
 
        1     STO 01
       14    STO 02
      314   STO 03    and if we choose   N1 = 10  &   N2 = 30

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

Notes:

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

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

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


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

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


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

   
      ( 241 bytes / SIZE 018 )

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

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

Example1:
     Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = -0.1
 
        1     STO 01
       14    STO 02
      314   STO 03  
     -0.1   STO 04  and if we choose   N1 = 10  &   N2 = 30

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

  and  R00 = L = -0.00003033995617  Gy-2


Example2:      Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1
 
        1     STO 01
       14    STO 02
      314   STO 03  
      0.1    STO 04  and if we choose   N1 = 10  &   N2 = 30

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

  and  R00 = L = -0.00003032063771  Gy-2

Notes:

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

-Register R17 may be replaced by synthetic register M.

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


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

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


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

   
      ( 197 bytes / SIZE 015 )

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

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

Example1:
     Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = -0.1
 
        1     STO 01
       14    STO 02
      314   STO 03  
     -0.1   STO 04  and if we choose   N1 = 1  &   N2 = 60

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

  and  R00 = L = -0.00003033995617  Gy-2


Example2:      Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1
 
        1     STO 01
       14    STO 02
      314   STO 03  
      0.1    STO 04  and if we choose   N1 = 15  &   N2 = 60

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

  and  R00 = L = -0.00003032063771  Gy-2


Note:

-The following variant employs Carlson elliptic integrals:


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

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


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

   
      ( 225 bytes / SIZE 017 )

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

 
Example1:
     Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = -0.1
 
        1     STO 11
       14    STO 12
      314   STO 13  
     -0.1   STO 14  

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

  and  R00 = L = -0.00003033995617  Gy-2


Example2:      Rmin = 1 Gly  ,  R0 = 14 Gly  ,   Rmax = 314 Gly   ,   R3 = 0.1
 
        1     STO 11
       14    STO 12
      314   STO 13  
      0.1    STO 14
   
               XEQ "ACU"   >>>>    T = 15.69777975 Gy                  
                                      RDN     P = 994.7216216 Gy
                                      RDN     H = 66.98887585 km/s/Mpc
                                      RDN     q  = -0.039995458

  and  R00 = L = -0.00003032063771  Gy-2

Notes:

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

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


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

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



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

   
      ( 302 bytes / SIZE 019 )

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

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

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

    71   STO 00

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

-If we choose  N1 = 20  &  N2 = 60

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

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


Notes:

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

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

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


4°)  Another Change of Variable


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

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

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

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

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


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

   
      ( 184 bytes / SIZE 015 )

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

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

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

    1    STO 01
   14   STO 02
  314  STO 03

-If we choose  N1 = 20  &  N2 = 60

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

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


Notes:

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

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


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

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


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

   
      ( 206 bytes / SIZE 017 )

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

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

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

    1    STO 01
   14   STO 02
  314  STO 03

-If we choose  N1 = 8  &  N2 = 24

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

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



References:

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

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