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

 3°)  More General Universes


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

 4°)  Another Change of Variable


Latest Update:   paragraph 3°)c)

-These programs solve Einstein's equations

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

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

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

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

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


 R

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

1°)  Radiation Only - Negative Pressure
 

       a) Time <-> Scale Factor


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

-Einstein's equations become:

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

-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  "AUR"  returns  T , P , H , q  &  cosmological constant L  is stored in R00


Data Registers:             R00 = L 

                                        R01 = R2min           R04-R05: temp
                                        R02 = R2            
                                        R03 = R2max
Flags: /
Subroutines: /

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

 
      ( 82 bytes / SIZE 008 )

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

 
Example:
     Rmin = 1    R0 = 14   Rmax = 314

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

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


2°)  Zero-Pressure & Negative Density of Matter


       a) Period of the Universe


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

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

-Multiplying by R'  gives:

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

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

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

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

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

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

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

-The following routine employs Gauss-Chebyshev formula:


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

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


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

 
      ( 91 bytes / SIZE 008 )

 
            STACK             INPUT          OUTPUT
                 X                 /                  P


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

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

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

         100  XEQ 10  >>>>    P = 978.7265794  GY  

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

 

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

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


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

 
      ( 72 bytes / SIZE 008 )

 
            STACK             INPUT          OUTPUT
                 X                 N                  P

     Where  N = number of subintervals for Gauss-Chebyshev formula

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

-If you choose  N = 120

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

 And  R00 = L = -0.00003033029693

Note:

-Exact value = 993.8651296  ( rounded to 7 decimals )


       b) Time -> Scale Factor


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




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

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


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

 
      ( 144 bytes / SIZE 011 )

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

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

-Calculate  R(16)

   -0.0000314   STO 00

            0           STO 01
          1.4          STO 02
            0           STO 03

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

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

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


       c) Scale Factor -> Time


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

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

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

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


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

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


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

 
      ( 176 bytes / SIZE 014 )

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

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

-Calculate  t(15)

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

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

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

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


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

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


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

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



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



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

Data Registers:              R00 = L        

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

-Line 150 is a three-byte GTO 02


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

   
      ( 266 bytes / SIZE 016 )

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


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

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

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


Notes:

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


3°)  More General Universes


       a) Period of the Universe


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

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

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


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

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


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

 
      ( 89 bytes / SIZE 010 )

 
            STACK             INPUT          OUTPUT
                 X                 N                  P

     Where  N = number of subintervals for Gauss-Chebyshev formula

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

-If you choose  N = 120

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

 And  R00 = L = -0.00003033995617


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

-If you choose  N = 120

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

 And  R00 = L = -0.00003032063771

Note;

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


       b) Age + Period of the Universe


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

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

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


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

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


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

   
      ( 248 bytes / SIZE 017 )

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

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

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

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

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

Notes:

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

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

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


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

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


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

   
      ( 241 bytes / SIZE 018 )

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

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

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

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

  and  R00 = L = -0.00003033995617  Gy-2


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

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

  and  R00 = L = -0.00003032063771  Gy-2

Notes:

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

-Register R17 may be replaced by synthetic register M.

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


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

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


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

   
      ( 197 bytes / SIZE 015 )

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

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

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

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

  and  R00 = L = -0.00003033995617  Gy-2


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

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

  and  R00 = L = -0.00003032063771  Gy-2


Note:

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

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


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


Data Registers:              R00 = L        

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


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

   
      ( 270 bytes / SIZE 016 )

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


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

    SCI 6

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

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


Notes:

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

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



Data Registers:              R00 = L        

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


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

   
      ( 274 bytes / SIZE 018+??? )

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


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

    SCI 7

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

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


Notes:

-The precision is controlled by the display setting.

-The following variant employs Carlson elliptic integrals:


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

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


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

   
      ( 225 bytes / SIZE 017 )

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

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

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

  and  R00 = L = -0.00003033995617  Gy-2


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

  and  R00 = L = -0.00003032063771  Gy-2

Notes:

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

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

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

  and with this formula, the program becomes slightly shorter.


Data Registers:              R00 = L        

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


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

   
      ( 194 bytes / SIZE 014 )

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


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

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

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


Notes:

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


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

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


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

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



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


       c) My Favorite Focal Program


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

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

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

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


Data Registers:              R00 = L        

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

-Line 153 is a three-byte GTO 02


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

   
      ( 308 bytes / SIZE 016 )

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


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

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

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


Notes:

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

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


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


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

   R13 = (Omega)mat       =  -0.07704871226
   R14 = (Omega) lambda =  -0.002137513639
   R15 = (Omega)rad       =  0.0004953126994


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 )