Gravity & Time Manual


 Overview
 

-This module deals with Cosmology in general relativity, Earth gravity, geodesic distances and different time scales ( UT , TDB , TCB ).
-Ellipsoidal models of the Earth are used, spheroidal and scalene ellipsoids.

-A Luni-Solar calendar may also be used... just for fun !
 
 
 
XROM  Function  Desciption
 16,00
 16,01
 16,02
 16,03
 16,04
 16,05
 16,06
 16,07
 16,08
 16,09
 16,10
 16,11
 16,12
 16,13
 16,14
 16,15
 16,16
 16,17
 16,18
 16,19
 16,20
 16,21
 16,22
 16,23
 16,24
 16,25
 16,26
 16,27
 16,28
 16,29
 16,30
 16,31
 16,32
 16,33
 16,34
 16,35
 16,36
 16,37
 16,38
 16,39
 16,40
 16,41
 16,42
 16,43
 16,44
 16,45
 16,46
 16,47
 16,48
 16,49
 16,50
 16,51
 16,52
 16,53
 16,54
 16,55
 16,56
 16,57
 16,58
 16,59
 16,60
 16,61
 16,62
 16,63
 -GRVITME 1B
  ABC
  AF
  "AUM"
  "AUR"
  "BCT"
  "BDT"
  "COSMO"
  "COSMO+"
  "CSM"
  "DT"
  "DTTM"
  "ETP"
  "EUCL"
  "FWD"
  "GAB"
  "GAB2"
  "GEO"
  "GEO+"
  "GR-LS"
  "GRV"
  "GRV+"
  "GRV3"
  "GRV3+"
  "GRV4"
  "GRV4+"
  HUBBLE
  "J0"
  J2
  "J2TM"
  "LS-GR"
  "SPOL"
  "SPOL+"
  "TGD"
  "TGD2"
  "TGD3"
  "TOL"
  "TT-UT"
-AUX FNS
  1/Z
  CDAY
  CBRT
  CMD
  CROSS
  ELLOX+
  DOT
  FMD
  GW
  GX
  JDAY
  NORM
  "P2"
  "P3"
  "P4"
  R-S
  RFZ
  RJZ
  S-R
  ST<>A
  HMT
  LEG
  #GR
  X/E3
  X^3
 Section Header
 3 semi-axes of Triaxial Ellipsoid model of the Earth
 Equatorial radius & 1 / flattening of the Earth
 Age of the Universe ( Radiation + Dust )
 Age of a Universe ( Radiation without Dust )
 Barycentric Coordinate Time
 Barycentric Dynamical Time
 Redshift -> Several Distances
 Idem + prompts for inputs
 Redshift -> Several Distances in an Empty Universe
 Number of Days since 2000/01/01 -> Date ( YMD )
 Idem with TIME module
 Einstein's Twin Paradox
 Redshift -> Distance in an Euclidean Universe
 Geodesic Distance: Forward Problem
 An Example of Subroutine to get the Speed of Light
 2nd Example for the Speed of Light
 Geocentric Coordinates
 Idem + Prompts
 Example of a Luni-Solar Calendar 
 Gravity on a Scalene Ellipsoid
 Idem + Prompts
 Gravity on a Spheroidal Earth
 Idem + Prompts
 Gravity from the Geopotential ( only a few terms )
 Stores the constants used by "GRV4"
 Hubble's Constant
 Number of Days since 2000/01/01 + Day of the week
 Idem for Gregorian dates ( without DOW )
 Idem with TIME Module
 Example of a Luni-Solar Calendar
 Speed of Light in General Relativity / coord. time
 Idem + Prompt
 Terrestrial Geodesic Distance ( Andoyer's 2nd order )
 Terrestrial Geodesic Distance ( Vincenty's formula )
 Geodesic Distance - Triaxial Ellipsoidal Earth
 Redshift -> Distance in Tolman Universes
 Terrestrial Time minus Universal Time
 Section Header
 Inverse of a complex number
 Julian Day -> Date
 Cube Root
 Ceil ( Y Z / X ) ( M-code routine used by LuniSol cal. )
 Cross Product ( X Y Z ) x ( M N O )
 Loxodromics on an Ellipoidal Earth
 Dot Product ( X Y Z ) . ( M N O )
 Floor ( Y Z / X ) ( M-code routine used by LuniSol cal. )
 M-Code for Gauss-Legendre 6-point formula
 Idem
 Date -> Julian Day
 Norm of a 3D-vector ( X Y Z )
 Quadratic Equation a.x^2+b.x+c = 0
 Cubic Equation a.x^3+b.x^2+c.x+d = 0
 Quartic Equationx x^4+a.x^3+b.x^2+c.x+d = 0
 Rectangular-Spherical Conversion
 Carlson Elliptic Integral 1st kind RF(x,y+iz,y-iz)
 Carlson Elliptic Integral 3rd kind RJ(x,y+iz,y-iz,t>0)
 Spherical-Rectangular Conversion
 Exchanges X-Y-Z <> M-N-O
 Hermite Polynomial
 Legendre Polynomial
 An M-Code Subroutine called by "GRV3"
 X divided by 1000
 Cube

 
 

Three Semi-Axes of a Scalene Ellipsoidal Model of the Earth
 

 ABC  places

  a = 6378.171274 km  in X-register
  b = 6378.101946 km  in Y-register
  c = 6356.751868 km  in Z-register

-The previous contents of registers X Y Z are overwritten.
-The last decimals are somewhat meaningless !
 

Equatorial Radius of the Earth & Inverse of the Earth Flattening
 

 AF  returns

   a  = 6378.13661  km  in X-register
  1/f =  298.256421   in Y-register

-The previous contents of X & Y are saved in Z & T
-These values seem the most recent estimations.
 

Gravity on a Scalene Ellipsoid
 

-"GRV" employs the formulae given in reference [6] and the semi-axes above.
 
 
      STACK        INPUTS      OUTPUTS
           Z         L ( ° ' " )          g(0)
           Y         B ( ° ' " )          g(0)
           X          h ( m )          g(h) 

 
Example1:     L =  77°03'56" W ,  B = 38°55'17"2 N ,   h = 67m                ( U.S. Naval Observatory , Washington D.C. )

 -77.0356    ENTER^
  38.55172   ENTER^
       67         XEQ "GRV"  >>>>   g(67) = 9.800516275  m/s2
                                          RDN    g(0)  = 9.800723034  m/s2
 

Example2:     L = 116°51'50"4 W ,    B = 33°21'22"4 N    h = 1706 m                ( Mount Palomar Observatory )

 -116.51504    ENTER^
   33.21224      ENTER^
      1706            R/S     >>>>   g(1709) = 9.790660012  m/s2
                                     RDN       g(0)    = 9.795923287  m/s2
 

Notes:

-The longitudes are positive East.
-Latitudes & longitudes are both supposed geodetic, but using a geocentric longitude will not change the result significantly.
-The results are very accurate provided the altitude remains relatively small
-Unfortunately, the Earth gravitational field is much more complex than that of a triaxial ellipsoid,
 so this program only gives slightly better results than the classical Somigliana's formula.

-XEQ "GRV+" to get PROMPTs for the inputs.
 

Gravity on a Spheroidal Earth
 

-The gravity g above an ellipsoid of revolution may be computed by rigorous formulae ( cf reference [7] ):
 
 
      STACK        INPUTS      OUTPUTS
           Z             /         B ( deg )
           Y         b ( ° ' " )          R ( m )
           X          h ( m )       g ( m/s2 )

   where   R = distance from the center of the Earth
    and      B = geocentric latitude

Example1:     b = 38°55'17"2     h = 23456 m

  38.55172   ENTER^
    23456      XEQ "GRV3"  >>>>   g = 9.728751603 m/s2
                                           RDN    R = 6393194.708 m
                                           RDN    B = 38°73415862
 

Example2:     b = 38°55'17"2     h = 12345678 m

  38.55172   ENTER^
 12345678   XEQ "GRV3"  >>>>   g = 1.078713339 m/s2
                                            RDN   R = 18715394.22 m
                                            RDN   B = 38°85746754

Note:

 XEQ "GRV3+" to get PROMPTs for the inputs
 

Gravity from the Geopotential ( a few terms )
 

-"GRV4" uses the Earth's gravitational potential U to compute the gravity ( spherical harmonics )
-Thousands - or at least hundreds - of terms are needed to get full accuracy.
-This routine employs only a few terms of the geopotential which must be stored in R19 thru R54 by "GRV4+" before executing "GRV4"
 
 
      STACK        INPUTS      OUTPUTS
           Z          B (deg)             /
           Y          L (deg)             /
           X          R  (m)        g ( m/s2)

  where  R , L , B  are the geocentric coordinates.

      R = distance from the center of the Earth
      L = longitude,
      B = geocentric latitude ( not geographic latitude )

-"GEO" may be used to calculate R & B from the altitude and the geographic latitude.

Example:      R = 6369806 m    L = -77.065556°    B = 38.733471°

  XEQ "GRV4+" to get the constants in the proper registers, then the HP-41 displays

    38.733471   ENTER^
  -77.065556    ENTER^
     6369806     R/S  >>>>   g = 9.800330901 m/s2            ---Execution time = 79s---

         gR = R01 = -9.800278379  m/s2
         gL = R02 =   0.000091758  m/s2
         gB = R03 = -0.032085221  m/s2

Note:

-The constants in R19 to R54 are unchanged by this program,
 so you can XEQ "GRV4" instead of "GRV4+" to get another g-value at another location
 

Hubble's Constant
 

 HUBBLE  returns in X-register Hubble's constant: 1 /  H0 = 1.377172143   E10  years

-Most of the decimals are of course doubtful
-They have been chosen to give H0 = 71 km / s / Mpc

-The relation betxeen  H0 = x km/s/Mpc  &  1 / H0 = y  109 years is    x y = 977.79222143
 

Age of the Universe ( Dust + Radiation )
 

-There are many subcases and the following program only works if the cosmological constant is positive

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

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

-These conditions are precisely satisfied by our Universe !
 
 
      STACK        INPUTS      OUTPUTS
           Z             ¶0            k
           Y             L0            R0
           X      (Omega)mat            t0

                  t0 = age of the Universe ( in years )
   Where    R0 = Radius of the Universe ( in light-years )
                  k =  +1 , 0 , -1  for Spherical , Euclidean , Hyperbolic spaces

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

   0.000049  ENTER^
      0.519     ENTER^
      0.044     XEQ "AUM"  >>>>  t0 = 1.5587728 E10 years
                                                     R0 = 2.0833961 E10 light-years
                                                     k  =  -1  ( hyperbolic space )

Example2:      With    (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,  (Omega)rad  = 0.000049

    ( These Omega-values are close to those suggested by WMAP = NASA's Wilkinson Microwave Anisotropy Probe )

      0.000049  ENTER^
        0.732     ENTER^
        0.271     XEQ "AUM"  >>>>  t0 = 1.3667886 E10 years
                                                      R0 = 2.4940750 E11 light-years
                                                       k  =  +1  ( spherical space )

Example3:      With    (Omega)mat = 0.2679 ,  (Omega)lambda  = 0.732 ,  (Omega)rad  = 0.0001

    .0001  ENTER^
     .732   ENTER^
    .2679  XEQ "AUM"  >>>>  t0  = 1.3693298 E10 years
                                                R0 = 1.3771721 E10 light-years          R0 = c/H0    when the Universe is Euclidean ( conventionally )
                                                 k  =  0  ( Euclidean space )

Note:

-"AUM" calls "P4" ( quartic equations ) and the Carlson elliptic integrals are computed with the M-Code routines RFZ & RJZ
 

Age of the Universe ( Radiation without Dust )
 

-These kinds of Universe are also called Tolman Universes.
-The calculations are much easier and do not require elliptic integrals

-Here,  (Omega)mat = 0
 
 
      STACK        INPUTS      OUTPUTS
           Z              /            k
           Y             ¶0            R0
           X             L0            t0

                 t0 = age of the Universe ( in years )
   Where    R0 = Radius of the Universe ( in light-years )
                  k =  +1 , 0 , -1  for Spherical , Euclidean , Hyperbolic spaces

Example:       ¶0 =  (Omega)rad = 49 10 -6  ,  L0 = (Omega)lambda = 0.61

      49 E-6  ENTER^
        0.61   XEQ "AUR"  >>>>   t0 = 1.8236300 E10   years
                                        RDN  R0 = 2.2053788 E10   light-years
                                        RDN   k  = -1  ( hyperbolic Universe )
 

Redshift >>> Distance of a Galaxy
 

-This problem involves the same kind of integrals as "AUM"
-So, quartic equation + Carlson elliptic integrals
-The same restrictions apply
 

Data Registers:   R00 thru R32

      R10 =  z + 1                  R13 =  (Omega)rad
      R11 =  (Omega)mat       R14  = (Omega)tot - 1
      R12 =  (Omega)lambda  R15 thru R32:  temp

 -when the program stops:

             R00 = k       ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
             R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )

Flags:  F00  set flag F00  if you don't want to re-calculate the age of the Universe t0
              F01 & F02 are used by "P4"
 
 
      STACK        INPUTS      OUTPUTS
           T            ¶0            VR
           Z             L0            DL
           Y      (Omega)mat            D0
           X              z            D
           L              /            t0

        D   =  light-travel time distance ( in light-years ) = light-travel time ( in years )
        D0 = comoving distance ( in light-years )
        DL = luminosity-distance ( in light-years )
        VR = recessional velocity ( speed of light = 1 )
        t0  = age of the Universe ( in years )

Example:      With    (Omega)mat = 0.044 ,  (Omega)lambda  = 0.519 ,  (Omega)rad  = 0.000049   &   z = 7

      0.000049  ENTER^
        0.519     ENTER^
        0.044     ENTER^
            7        XEQ "COSMO"    >>>>   D  =  1.4119088  E10  l-y     = R01
                                                    RDN    D0 = 3.4190016  E10  l-y     = R02
                                                    RDN    DL = 4.1392272  E11  l-y     = R03
                                                    RDN    VR = 2.482624772               = R04         ( unit c = 1 )
                                                  LASTX   t0  = 1.5587728  E10  years = R05

   and   R00 = k = -1  ( hyperbolic space )
           R06 =  R0 = 2.0833961  E10  light-years

Notes:

-XEQ "COSMO+" to get PROMPTs for the 4 inputs
-Execution time is about 1 minute.
 

Redshift >>> Distance of a Galaxy in an Empty Universe
 

-The problem is of course much simpler in an empty universe.
-"CSM" only works with a cosmological constant such that  0 < (Omega)lambda  < 1
 
 
      STACK        INPUTS      OUTPUTS
           T             /            VR
           Z             /            DL
           Y             L0            D0
           X              z            D
           L              /            t0

        D   =  light-travel time distance ( in light-years ) = light-travel time ( in years )
        D0 = comoving distance ( in light-years )
        DL = luminosity-distance ( in light-years )
        VR = recessional velocity ( speed of light = 1 )
        t0  = age of the Universe ( in years )

Example:       (Omega)lambda  = 0.519    &   z = 7

        0.519     ENTER^
            7        XEQ "CSM"    >>>>   D  =  1.4892172  E10  l-y
                                              RDN    D0 = 3.7410987  E10  l-y
                                              RDN    DL = 5.1055493  E11  l-y
                                              RDN    VR = 2.716507697                   ( unit c = 1 )
                                            LASTX   t0  = 1.7367386  E10  years

Note:

 k = -1  for these models
 

Redshift >>> Distance of a Galaxy in a Euclidean Universe
 

-Given the cosmological parameter and the redshift z
  "EUCL" returns the light-travel time distance D and the age of the Universe.
 
 
      STACK        INPUTS      OUTPUTS
           Y             L0            t0
           X              z            D

    D   =  light-travel time distance ( in light-years ) = light-travel time ( in years )
     t0  = age of the Universe ( in years )

Example:     L0 = 0.732  &  z = 7

    0.732  ENTER^
       7      XEQ "EUCL"   >>>>   D = 1.2915891  E10 ly
                                       X<>Y   t0 = 1.3698976  E10 years
 

Redshift >>> Distance of a Galaxy in Tolman Universes
 

-Like with"AUR",  (Omega)mat = 0
-Such models have been studied by Tolman.

-There are many cases, and "TOL" only works for Big Bang Universes with a positive Cosmological Constant.
-The comoving radial distance involves elliptic integrals and "TOL" calculates the light-travel time Distance D,
  the age of the Universe t0 , the current radius of the Universe R0 and the constant  k
  ( D & t0 may be expressed in terms of elementary functions )
 
 
      STACK        INPUTS      OUTPUTS
           T              /            k
           Z             ¶0            R0
           Y             L0            t0
           X              z            D

Example:       ¶0 =  (Omega)rad = 49 10 -6  ,  L0 = (Omega)lambda = 0.61  ,   z  = 7

      49 E-6  ENTER^
        0.61    ENTER^
           7      XEQ "TOL"  >>>>  D = 1.5726496  1010   light-years
                                        RDN   t0 = 1.8236300  1010   years
                                        RDN  R0 = 2.2053788 1010   light-years
                                        RDN   k  = -1  ( hyperbolic Universe )
 

Barycentric Dynamical Time
 

-This program calculates TDB-TT. It employs the first terms of the series given in reference [8],
  those whose amplitude > 0.4 second over the time-span [1000,3000]
-The errors should not exceed 1 or 2 microseconds over this interval of time.
 
 
      STACK        INPUTS      OUTPUTS
           Y  YYYY.MNDD             /
           X      HH.MNSS     TDB-TT (s)

 
Example1:      2012/08/29  16h41m37s  ( TT or TDB )

    2012.0829    ENTER^
        16.4137    XEQ "BDT"  >>>>   TDB-TT =  - 0.0013227 s
 

Example2:      3000/04/04  1h23m45s  ( TT or TDB )

    3000.0404   ENTER^
          1.2345   XEQ "BDT" ( or  RTN  R/S )    >>>>    TDB-TT = + 0.0015386 s
 

Barycentric Coordinate Time
 

-Anothe time scale may also be used:  TCB = Temps Coordonné Barycentrique = Barycentric Coordinate Time.
-Whereas  TDB - TT  remains smaller than 0.002 second at least several millenia aroud J2000, the difference  TCB - TDB  is much larger:
-It was almost 0 in 1977 but in 2012 , TCB -TDB is about 17 seconds.

Formula:

   TCB = TDB + 1.550519768 E-8 ( JDTCB - 2443144.5003725 ) 86400 + 6.55 E-5 seconds

     where  JDTCB  is the Julian Date in the TCB scale.

-"BCT" uses an iteration to compute  TCB - TDB  for a given date expressed in  TDB
 
 
      STACK        INPUTS      OUTPUTS
           Y  YYYY.MNDD             /
           X      HH.MNSS    TCB-TDB (s)

  Where the input time is expressed in TDB

Example:       3000/04/04  1h23m45s  TDB

     3000.0404   ENTER^
           1.2345   XEQ "TCB"   >>>>   500.6752393 seconds
 

Terrestrial Time minus Universal Time
 

-This routine uses the polynomials given in reference [9]
 
 
      STACK        INPUTS      OUTPUTS
           X          year     TT-UT (sec)

 
Examples:

 -2000   XEQ "TT-UT"  >>>>  46675.68 s
   400    XEQ "TT-UT"  >>>>   6699.22 s
  1200   XEQ "TT-UT"  >>>>    736.44 s
  1680   XEQ "TT-UT"  >>>>     15.31 s
  1760   XEQ "TT-UT"  >>>>     14.87 s
  1841   XEQ "TT-UT"  >>>>      5.52 s
  1880   XEQ "TT-UT"  >>>>    -5.01 s
  1906   XEQ "TT-UT"  >>>>      5.10 s
  1934   XEQ "TT-UT"  >>>>     23.86 s
  1951   XEQ "TT-UT"  >>>>     29.47 s
  1984   XEQ "TT-UT"  >>>>     53.73 s
  2000   XEQ "TT-UT"  >>>>     63.86 s
  2041   XEQ "TT-UT"  >>>>     85.52 s
  2100   XEQ "TT-UT"  >>>>    202.74 s
  3000   XEQ "TT-UT"  >>>>   4435.68 s

Note:

-In most cases, the program does not stop at the END
-So, press XEQ "TT-UT" or  RTN  R/S  for another date
 

Date <> Number of Days since 2000/01/01
 

-"J0" computes the number of days since 2000/01/01  0h ( i-e the Julian Day number minus 2451544.5 ) in X-register and the day of the week in T-register.
-Registers Y and Z are saved.
-The Gregorian calendar reform is taken into account: The day following 1582/10/04 is 1582/10/15
-"J0" ( and "DT" = the reverse operation ) are valid at least since 4713 B.C. up to ... the next calendar reform!
 
 
        STACK          INPUTS       OUTPUTS
             T               /           dow
             Z               Z             Z
             Y               Y             Y
             X   YYYY.MNDDdd             N
             L               /   YYYY.MNDDdd

  ( dow = day of week = 0 for Sunday, 1 for Monday, ... , 6 for Saturday )

-N= the number of days between a date YYYY.MMDDdd and 2000/01/01 at 0h  = the Julian Day number minus 2451544.5

Examples:

  April 4th 2134 at 6h AM    2134.040425  XEQ "J0"  >>>     N = 49036.25   RDN  RDN   dow = 0 = Sunday

    1234.0428    XEQ "J0"   >>>   N = -279651  RDN  RDN  dow = 5 = Friday
   -4123.0707   XEQ "J0"   >>>   N = -2236225  RDN  RDN  dow = 1 = Monday
 

-"DT" performs the reverse transformation
 
 
        STACK          INPUTS       OUTPUTS
             Z              Z             Z
             Y               Y             Y
             X               N   YYYY.MNDDdd
             L               /             N

 
Examples:

   N =   49036.25    XEQ "DT"  >>>   2134.040425
   N = -279651        XEQ "DT"  >>>   1234.0428
   N = -2236225      XEQ "DT"  >>>  -4123.0707

Notes:

 J2 is an M-Code routine that takes a Gregorian date ( even before 1582 ) YYYY.MNDD and returns the number of days since 2000/01/01
 J2TM  does the same thing but with the TIME module functions

 Like "DT",  DTTM  takes the number of days since 2000/01/01 and returns the Gregorian date YYYY.MNDD
 

Einstein's Twin Paradox
 

-"ETP" uses Einstein's special relativity to deal with uniformly accelerated linear motion of a spacecraft.
-In special relativity, the speed is also  v = dx/dt , but acceleration is  a = ( dv/dt ) ( 1 - v2/c2 )-3/2   where  t = time passed on Earth
-The proper Time T ( passed on board ship ) verifies  dT = dt ( 1 - v2/c2 )1/2

-We use the following units:

    -Times ( t , T ) are expressed in Julian years ( 365.25 days )
    -Distance  D is expressed in light-years ( one light-year = 9.460,730,472,580,8 1015 m. )
    -The unit of accelerations ( a ) is the "standard" gravity on Earth  g  = 9.80665 m/s2

    -The velocity of light = c = 299,792,458 m/s  ( from the definition of 1 meter )
 

*********************************************************************************************

                                                    Return Journey
              <<<<<Deceleration "a" <<<<<  |  <<<<<Acceleration "-a" <<<<<<<
 Earth=O-------------------------------D/2------------------------------------D=(Star)--------------->   x
              >>>>>Acceleration "a" >>>>>  |  >>>>>Deceleration "-a" >>>>>>>
                                                  Outward Journey

-The graphic representation of  x(t) looks like this:

        x
         |
         | D                                              *
         |                                    *                         *
         |                             *                                       *
         | D/2                 *                                                  *
         |                  *                                                              *
         |           *                                                                            *
         |*---------------|----------------|-----------------|-----------------*-----  t
    Takeoff                          Spacecraft reaches the star                        Landing
 

*************************************************************************************************
 
 
      STACK        INPUTS      OUTPUTS
           T             /        1-vmax/c
           Z             /          vmax/c
           Y             /             tE
           X             /            TS

       TS = time passed on board ship
        tE = time passed on Earth.
      vmax = maximum speed  ( when x = D/2 )

Example1:    On 2007/12/25  immediate boarding for the first interstellar flight in the direction of Toliman ( alpha Centauri ) Distance D = 4.343 light-years.
                           Fox Mulder remains on Earth but Dana Scully's flying saucer takes off at 10h13m TT , acceleration = 0.7g

   XEQ "ETP"   >>>>

 4.343  ENTER^
   0.7       R/S             >>>>      TS    =   8.837390060 years
                                 RDN        tE    = 13.09998313  years
                                 RDN    vmax/c  =  0.9211383858              vmax  = 92.11383858% * 299792458 m/s = 276150341 m/s
                                 RDN  1-vmax/c = 0.07886161418

-When she steps from the ship, Dana's watch shows 2016/10/26 , 6h46m41s
  wheras according to Mulder, it's 4h40m08s o'clock and we are on 2021/01/30
 

Example2:     With a = 1g and D = 2,200,000 light-years   ( approximately the distance of the Andromeda galaxy )   "ETP" returns:

         TS    =      56.71150062    years
          tE    =     4,400,003.874  years
      vmax/c =        1
   1-vmax/c = 3.877715920 E-13  whence  vmax/c  was actually  0.999,999,999,999,612,228,408

Note:

-Caculating   1-vmax/c is only useful if the maximum  v/c  is very close to 1
 

Speed of Light in General Relativity / coordinate time
 

-If the clocks are synchronized along the path of the light and if we use the proper time, the speed of light is constant.
-In this case, no program is needed to calculate it:  c = 299792458 m/s

-The following routine employs the coordinate-time t to measure the velocity of light  V = dL/dt
-We may have V < c or V > c  or V = c

-The propagation of light is characterized by null geodesics:

        ds2 = gab dxa dxb = 0       a , b  =  0 , 1 , 2 , 3    where    x0 = c.t  &  x1 , x2 , x3  are spatial coordinates

-This may be re-written

              ds2 = g00 dx0 dx0 + 2 g0i dx0 dxi + gij dxi dxj       i , j = 1 , 2 , 3
                    = [ sqrt(g00) dx0 + g0i dxi / sqrt(g00) ]2 -  ( g0i  g0j / g00 ) dxi dxj  + gij dxi dxj
                    =  g00 ( dx0 )2 [ 1 + ( g0i / g00) ( dxi / dx0 ) ]2  -  dL2

   where    dL2 =  hij dxi dxj    is the spatial metric  with  hij = - gij + ( g0i  g0j / g00 )

 [ The spatial tensor is usually denoted  gij  ( the geek letter gamma ), but I've prefered to use hij  in case your browser ignores greek symbols ]

-The coordinates of the velocity of light are  Vi = dxi / dt  and its modulus is  V = dL/dt = sqrt ( hij Vi Vj )

   whence  ds2 / dt2 = 0 =  g00 c2 [ 1 + ( g0i / g00) ( 1/c ) ( dxi / dt ) ]2  -  V2  =  g00 c2 [ 1 + ( g0i / g00) Vi / c ]2  -  V2
 

-The direction of propagation is defined by a 3-vector  ki = dxi / dL = ( dxi / dt ) ( dt / dL ) = Vi / V      (  hij ki kj = 1 )

  so,   V = c Sqrt(g00)  +  ( g0i / sqrt(g00) ) ki V

-Finally, it yields:        V/c  = sqrt (g00) / [ 1 - ki g0i / sqrt(g00) ]
 

Data Registers:           •  R00 = function name                        ( Registers R00 thru R07 are to be initialized before executing "SPOL" )

                                      •  R01 = x1                   •  R05 = k1           R08 = g00       R12 = g01       R15 = g12       R17 = g23
                                     •  R02 = x2                   •  R06 = k2           R09 = g11       R13 = g02       R16 = g13
                                      •  R03 = x3                   •  R07 = k3           R10 = g22       R14 = g03
                                      •  R04 = x0 = c.t                                        R11 = g33

                                                         R09 = h11       R13 = h02       R16 = h13
                                                         R10 = h22       R14 = h03
                                                         R11 = h33
Flags: /
Subroutine:   A program that takes  x1 ,  x2  ,  x3  ,  x0  in registers  R01-R02-R03-R04  and calculates
                       and stores the 10 components gab  into  R08 to R17  as shown above.
 
 
      STACK        INPUT      OUTPUT
           X            /           v/c

   where  v  is the speed of light measured with the time coordinate

Example1:          ds2 = [ 1 - 1 / (4r) + r2 / 1000 ] c2 dt2 - ( dx2 + dy2 + dz2 ) / [ 1 - 1/(4r) + r2/1000 ] - ( 2 / r3 ) ( y.dx - x.dy ) c dt

   with  r = sqrt ( x2 + y2 + z2 )

  •   Evaluate the speed of light at a point  P(1,2,3,0)

          a)  in the direction defined by the vector k(4,5,6)
          b)  in the direction defined by the vector k(0,1,0)

-We have  g00 = 1 - 1/(4r) + r2/1000 = - 1 / g11 =  - 1 / g22 =  - 1 / g33   ,  g01 = - y / r3  ,  g02 = x / r3  , the other components of the metric tensor = 0
 

>>>  The subroutine "GAB" is already in the module and stores the components of this metric tensor in the proper registers
 

  "GAB"   STO 00                                     ( use XRQ "SPOL+" and you will get PROMPTs for all the inputs )

      1        STO 01
      2        STO 02
      3        STO 03
      0        STO 04

 a)  In the 1st direction

      4        STO 05
      5        STO 06
      6        STO 07

   XEQ "SPOL"  >>>>   V/c = 0.966923596

  b) In the 2nd direction

     0        STO 05
     1        STO 06
     0        STO 07

   XEQ "SPOL"  or  XEQ 10  ( faster )   >>>>  V/c = 0.992171327
 

  •   Evaluate the speed of light at a point  P(7,8,9,0)  in the direction defined by the vector k(2,3,4)

      7  STO 01      2  STO 05
      8  STO 02      3  STO 06
      9  STO 03      4  STO 07
      0  STO 04

  XEQ "SPOL"  >>>>   V/c = 1.084831634

Notes:

-XEQ "SPOL+" to get PROMPTs for the inputs

-The components of k may be multiplied by a positive constant without changing the results.

-Line 48  LBL 10  is only useful to get faster results at the same point.

-The metric tensor of this example is not very realistic though it could represent approximately the gravitational field of a rotating sphere
  in a universe with a ( very ) large negative cosmological constant.
 

Example2:          ds2 = [ 1 -  r2 / 10000 ] c2 dt2 - ( dx2 + dy2 + dz2 )  - ( 0.02  ) ( x.dy - y.dx ) c dt

   with  r = sqrt ( x2 + y2 + z2 )

  •   Evaluate the speed of light at the point  P(1,0,0,0)  in the directions  k(0,1,0)  and k(0,-1,0)

-We have  g00 = 1 -  r2/10000  ,   g11 = g22 = g33 = -1  ,  g01 = 0.01 y  ,  g02 = -0.01 x  , the other components of the metric tensor = 0
 

>>>  The subroutine "GAB2" is already in the module and stores the components of this metric tensor in the proper registers
 

 "GAB2"   STO 00                                     ( use XRQ "SPOL+" and you will get PROMPTs for all the inputs )

      1        STO 01     0   STO 05
      0        STO 02     1   STO 06
      0        STO 03     0   STO 07
      0        STO 04

  XEQ "SPOL"  >>>>   V/c = 0.990049504

     1  CHS  STO 06   R/S  or  XEQ 10  >>>>  V'/c = 1.010050504
 

Notes:

-This metric corresponds to a rotating frame of reference.
-The difference between V and V' explains the Sagnac effect ( interference )

-More generally, if the angular velocity = w , the metric may be written in cylindrical coordinates:

  ds2 = ( 1 - r2w2/c2 )  c2 dt2 - dr2 - r2 df2 - dz2 - 2 w r2 df dt

-And the spatial metric is

  dL2 = dr2 + dz2 + r2 df2 / ( 1 - r2w2/c2 )

-The ratio  circumference / radius = 2.p / ( 1 - r2w2/c2 ) 1/2 > 2.p

-As for the speed of light ( if its path is the circumference )   V/c = 1 ± r.w/c    ( if we neglect the terms of order 2 and higher )

-Note that we've just gotten a more accurate result in the example above with  r.w/c = 0.01
 

Remarks:

-"SPOL" is essentially useful if there is at least an index i  for which  g0i # 0
-Otherwise, the formula is much simpler:

     V/c = Sqrt(g00)

-In this case, V does not depend on the direction k and calculating the other  gab  is unuseful.
-For instance, with the Schwarzschild metric:

     ds2 = [ 1 - 2 GM / ( c2 r ) - L r2 / 3 ] c2 dt2 - ( dx2 + dy2 + dz2 ) / [ 1 - 2 GM / ( c2 r ) - L r2 / 3 ]    where   L  is the cosmological constant,

     G = gravitational constant = 6.673 E-11  m3/ kg / s2  and   M = mass of the central body ( for example, the Sun  M = 1.989 E30 kg )

-The speed of light is     V/c = [ 1 - 2 GM / ( c2 r ) - L r2 / 3 ] 1/2

-Neglecting the cosmological constant, the speed of light near the Sun ( at a distance = Sun's radius = 6.96 E8 m ) is

    V/c ~  0.999997878   whence   V ~  299791822 m / s
 

Example of a Luni-Solar Calendar
 

 "GR-LS" & "LS-GR" perform the conversion between Gregorian dates and a hypothetic Luni-Solar calendar
  that employs the following approximations:

   1 Lunar month = 334995 / 11344 days ~ 29.53058886  days
   1 tropical year / 1 lunar month = 774439 / 62615 ~ 12.36826639

which leads to  1 tropical year ~ 365.2421896 days
 
 
           STACK            INPUT1          OUTPUT1           INPUT2         OUTPUT2
               X     YYYY.MNDDGr         yyyy.mnddLS         yyyy.mnddLS    YYYY.MNDDGr
 
 
Examples:

  -3101.0123   XEQ "GR-LS"   >>>>      0.0101      XEQ "LS-GR"  ( or simply  R/S )  >>>> -3101.0123
   1979.0716   XEQ "GR-LS"   >>>>   5080.0721   XEQ "LS-GR"  ( or simply  R/S )  >>>>   1979.0716
   2012.1221   XEQ "GR-LS"   >>>>   5113.1308   XEQ "LS-GR"  ( or simply  R/S )  >>>>   2012.1221
   2015.1014   XEQ "GR-LS"   >>>>   5116.1001   XEQ "LS-GR"  ( or simply  R/S )  >>>>   2015.1014

Notes:

-Here, the years yyyy may also be negative ( otherwise, delete lines 48-40-39-38 ).
-These programs take the onset of the Kali-Yuga as origin i-e  -3101.0123 ( Gregorian Calendar )  -3101 = 3102 B.C.
 

Loxodromics on an Ellipoidal Earth
 

Formulae:

    L2 - L1 =  ( Tan µ ).{ Ln [ ( Tan ( 45° + b2/2 ) ) / (Tan  ( 45° + b1/2 ) ) ] - (e/2).Ln [ ( 1 + e.sin b2 ).( 1 - e.sin b1 )/( 1 - e.sin b2 )/( 1 + e.sin b1 ) ] }
            D =  [ a.( 1 - e2 )/ cos µ ] §b1b2  ( 1 - e2 sin2 t ) -3/2 dt     ( § = Integral )

-One could use any integration formula but "ELLOX" evaluates this integral by a series expansion up to the terms in  e6
  ( the first neglected term is proportional to  e8 )

            D = [ a.( 1 - e2 )/ cos µ ] [ ( 1 + 3e2/4 + 45e4/64 + 175e6/256 ).( b2 - b1 ) - ( 3e2/8 + 15e4/32 + 525e6/1024 ).( sin 2b2 - sin 2b1 )
                   + ( 15e4/256 + 105e6/1024 ).( sin 4b2 - sin 4b1 ) - ( 35e6/3072 ).( sin 6b2 - sin 6b1 ) ] + .......................

-However, we cannot use this formula if   cos µ = 0 , in this case ( i-e  if  b1 = b2 ) we have:

           D = a.( cos b ).( L2 - L1 ) / ( 1 - e2 sin2 b ) -1/2    the loxodromic line is the parallel of latitude b = b1 = b2
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )        µ ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

 
Example1:

     XEQ "ELLOX+"                L1^B.1=?

   -77.03560  ENTER^
    38.55172     R/S                 L2^B.2=?

      2.20138  ENTER^
    48.50112     R/S        >>>>   D = 6453.389603 km
                                     X<>Y    µ = 80°10'15"31
 

-Compare this value with the length of the corresponding geodesic:  6181.621427 km
-Rhumb-sailing is easier but slower than orthodromy.
 

Example2:

     XEQ "ELLOX+"                L1^B.1=?

           0   ENTER^
          30     R/S                      L2^B.2=?

        120  ENTER^
         30     R/S        >>>>   D = 11578.35295 km
                               X<>Y    µ = 90°
 

Notes:

-Due to roundoff errors, there is a loss of significant digits if the 2 latitudes are very close but not equal.
-Otherwise, the accuracy is of the order of a few centimeters.

-In order to avoid a DATA ERROR  & an OUT OF RANGE if the latitudes are +/- 90°
 key in  89°59'59"9999 instead of 90° , the difference is negligible.
 

Terrestrial Geodesic Distance ( Andoyer's 2nd order )
 

-Andoyer's formula gives a good accuracy.
-Moreover, the method in non-iterative and therefore, fast !
-It employs a spheroidal model of the Erath ( Ellipsoid of revolution )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )             /
           X       b2 ( ° ' " )        D ( km )

 
Example:

   77.0356    ENTER^
   38.55172  ENTER^
   -2.20138   ENTER^
   48.50112  XEQ "TGD"  >>>>   D = 6181.621449 km

-The error smaller than 15 mm !
-The accuracy is very good ( less than 1 meter ), provided D is not too large
-So, avoid nearly antipodal points...
 

Terrestrial Geodesic Distance ( Vincenty's formula )
 

-Vincenty's formula gives even more accurate results but, on the other hand, it's an iterative method.
-The forward and backward azimuths are also calculated.

-The longitudes are positive East
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )     Az2 ( ° ' " )
           Z       b1 ( ° ' " )     Az2 ( ° ' " )
           Y       L2 ( ° ' " )     Az1 ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

 where   Az1  =  forward azimuth  ( from P1 to P2 )
   and     Az2  =  back azimuth  ( from P2 to P1 )

Example:

   -77.0356    ENTER^
    38.55172   ENTER^
      2.20138   ENTER^
   48.50112  XEQ "TGD2"  >>>>   the successive Ln-values are displayed  and  finally

                                           D   =  6181.621427 km
                             RDN    Az1 =  51°47'36"8132
                             RDN    Az2 = -68°09'58"9654

-We also have  R06 = µ = azimuth at the equator = 37.74571892°

Notes:

-The accuracy is excellent, provided the 2 points are not nearly antipodal.
-Otherwise, the algorithm does not converge ! It happens when | Ln | exceeds 180°
 

Geodesic Distance: Forward Problem
 

-Given the coordinates ( L1 , b1 ) of the first point  P1 , the forward azimuth Az1 and the length D of the geodesic,
 "FWD" computes the longitude and the latitude ( L2 , b2 ) of the second point  P2
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )       b2 ( ° ' " )
           Z       b1 ( ° ' " )       b2 ( ° ' " )
           Y     Az1 ( ° ' " )       b2 ( ° ' " )
           X        D ( km )       L2 ( ° ' " )

 where   Az1  is the forward azimuth in the direction from the first point P1 towards the second point P2

Example:     L1 = 10°30'     b1 = 49°41'     Az1 = 12°24'     D = 16000 km       Calculate  L2 & b2

   10.30   ENTER^
   49.41   ENTER^
   12.24   ENTER^
  16000   XEQ "FWD"  >>>  the successive sigma-approximations are displayed and eventually, it yields:

               L2 = -177°03'08"003     X<>Y    b2 =  -14°06'40"7483

-This method works for antipodal points too and - more generally - for any length D !
 

Geodesic Distance - Triaxial Ellipsoidal Earth
 

 "TGD3" works for a triaxial ellipsoid
-The geodesic distance involves integrals which are approximated by Gauss-Legendre 6-point formula
-The method calculates 1 ( SF 01 ) or 2 terms ( CF 01 ) of Fourier series.

-Flags F03 & F04 are also used as shown below:
 

Flags:  F01-F03-F04

   CF 01  2 terms in the Fourier series
   SF 01  only 1 term is found

   SF 04   geocentric longitudes and latitudes

      CF 04 & CF 03  geocentric longitudes &  geodetic latitudes
      CF 04 & SF 03   geodetic longitudes and latitudes
 
 
      STACK        INPUTS      OUTPUTS
           T       l1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       l2 ( ° ' " )       D1 ( km )
           X       b2 ( ° ' " )  D2 or D1 ( km )
           L             /       D0 ( km )

          CF 01  ---Execution time = 3mn41s---
          SF 01   ---Execution time = 2mn14s---

  D0 is computed without Fourier series
  D1 is computed with 1 term in the Fourier series  ( SF 01 )
  D2 is computed with 2 terms in the Fourier series  ( CF 01 )

Example1:       Geodesic distance between the U.S. Naval Observatory at Washington (D.C.) and the "Observatoire de Paris"

   l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
   l2 = 2°20'13"8 E  ,    b2 = 48°50'11"2 N

    CF 01  CF 03  CF 04

  -77.0356    ENTER^
   38.55172  ENTER^
    2.20138   ENTER^
   48.50112  XEQ "TGD3"  >>>>   D2 = 6181.625506 km = R01
                                           X<>Y  D1 = 6181.625507 km
                                        LASTX   D0 = 6181.628094 km = R21

>>> Compared to the value given by the rigorous Jacobi's method ( 6181.625476 ),  D2 is in error by only 3 centimeters !
 

Example2:      Geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

   l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
   l2 = 18°28'35"7 E  ,  b2 = 33°56'02"5 S

    CF 01  CF 03  CF 04

  -77.0356    ENTER^
   38.55172   ENTER^
   18.28357   ENTER^
  -33.56025  XEQ "TGD3"  >>>>    D2 = 12709.56468 km = R01
                                           X<>Y   D1 = 12709.56609 km
                                        LASTX   D0 = 12709.56636 km = R21
 

Example3:    With ( -40° , -40° ) ( +120° , +50° ) it yields:

  •  If  CF 01  CF 03  CF 04      geocentric longitudes, geodetic latitudes

                  D2 = 18095.48412 km = R01
     X<>Y   D1 = 18095.48948 km
   LASTX   D0 = 18095.54448 km = R21

  •  If  CF 01  SF 03  CF 04      geodetic longitudes and latitudes

                  D2 = 18095.49331 km = R01
     X<>Y   D1 = 18095.49868 km
   LASTX   D0 = 18095.55361 km = R21

  •  If  CF 01  SF 04                  geocentric longitudes and latitudes

                  D2 = 18099.69052 km = R01
     X<>Y   D1 = 18099.69587 km
   LASTX   D0 = 18099.75037 km = R21

Notes:

-Longitudes are positive East.

-The accuracy of D0 is about 60 meters provided the points are not nearly antipodal.
-If flag F01 is set, D2 is not computed ( faster but less accurate in general ).
-The results should verify:  D2 <= D1 <= D0   and   D2  should be the most accurate result.
-However, the formulae don't make the differences between maximum and minimum: they just return an approximate extremum with  D1 & D2

-Unfortunately, the surface of the Earth is much more comlex than that of a scalene ellipsoid.
 

Geocentric Coordinates
 

-If the latitude b and the height h of an observer O are known, "GEO" calculates
 the geocentric latitude b' and the distance rho = OC  where C is the center of the Earth.
 

Formulae:

      (rho)  sin b' = a(1-f) sin u + h sin b                  a = 6378.13661 km
      (rho) cos b' =  a cos u + h cos b                      f = 1/298.256421

      where   tan u = (1-f) tan b
 
 
      STACK        INPUTS      OUTPUTS
           Y         h (km)        rho (km)
           X        b ( ° ' " )        b' ( ° ' " )

 
Example:   Find the geocentric coordinates of the Mount Palomar Observatory     b = 33°21'22"4     h = 1706 m = 1.706 km

   1.706      ENTER^
 33.21224  XEQ "GEOC"  >>>>    b'  = 33°10'47"124
                                         X<>Y   rho = 6373.415625 km

Note:

 XEQ "GEO+" to get PROMPTs for the inputs ... and the outputs !
 

Date <> Julian Day
 

JDAY and CDAY are mutually reciprocal date functions to convert a given calendar date into the Julian day number and back to the calendar date. Use flag 00 to select either Julian or Gregorian calendars in the conversions. The date format is also MM.DDYYYY regardless of the time module settings if there’s one.
 
 
       STACK  INPUT/OUTPUT  OUTPUT/INPUT
            X    MN.DDYYYY     Julian Day Nb

 
Example:  the date May 21st, 2014 corresponds to 2,456,799 (Gregorian calendar) or 2,456,812 (Julian calendar) day numbers.

You can also use JDAY to calculate the elapsed number of days between two dates – simply converting both to their Julian day numbers and subtracting them. If you do that you’ll notice a small discrepancy (18 days) between this approach and the results from DAYS – which leads me to believe that DAYS has some different convention – perhaps it assumes a 30-day month for financial calculations -  but unfortunately it appears to be a stealth function, as there is no documentation for it at all in the Securities Pac manual.

These two functions are based on the PPC routines JC and CJ – ported to an all-MCODE implementation to make it faster.
The formulas are given in the PPC ROM manual.
 

Rectangular<>Spherical Conversion
 

       x = r cos b cos l
       y = r cos b sin l
       z = r sin b

   where    x , y , z = rectangular coordinates,    r  = ( x2 + y2 + z2 )1/2  ,   l  = longitude ,  b = latitude
 
 
      STACK       INPUTS     OUTPUTS
           T            T            T
           Z            z            b
           Y            y            l
           X            x            r
           L            /      (x2+y2)1/2

 
Example:     x = 3 ; y = 4 ; z = -7   find the spherical coordinates ( in DEG mode )

    -7  ENTER^
     4  ENTER^
     3  XEQ "R-S"    r  =  8.602325267
                  RDN     = 53.13010235°
                  RDN     b = -54.46232221°

-The reciprocal function work similarly
 
 
      STACK       INPUTS     OUTPUTS
           T            T             T
           Z            b             z
           Y            l             y
           X            r             x
           L            /      (x2+y2)1/2

 
Example:     r = 10 ; l = 124° ; b = 37°   find    x ; y ; z   ( in DEG mode )

   37   ENTER^
  124  ENTER^
   10  XEQ "S-R"   x = -4.465913097
                  RDN   y =  6.620988446
                  RDN   z =  6.018150232

Note:

 R-S and S-R work in all angular modes
 

Hermite Polynomial
 

Formulae:      Hn(x) = 2x.Hn-1(x) - 2(n-1).Hn-2(x)  ;  H0(x) = 1  ;  H1(x) = 2x
 
 
      STACK      INPUTS      OUTPUTS
           Y            n             n
           X            x          Hn(x)
           L            /             x

 
Example:    Calculate   H7(4.9)

    7       ENTER^
  3.14  XEQ "HMT"  >>>>  H7 (3.14) =   73726.24332
 

Legendre Polynomial
 

Formulae:      n.Pn(x) = (2n-1).x.Pn-1(x) - (n-1).Pn-2(x)  ;  P0(x) = 1  ;  P1(x) = x
 
 
      STACK      INPUTS      OUTPUTS
           Y            n             n
           X            x          Pn(x)
           L            /             x

 
Example:    Calculate   P7(4.9)

    7    ENTER^
  4.9   XEQ "LEG"  gives   P7(4.9) =  1698444.018
 

Quadratic Equation
 

"P2" solves the equation:   a.x2+b.x+c = 0    with a # 0
 
 
 
 
      STACK        INPUTS      OUTPUTS
           Z          a # 0           c/a
           Y             b             y
           X             c             x

-If  CF 00  the 2 solutions are  x , y
-If  SF 00   ------------------  x+i.y , x-i.y

Examples:     Solve  2.x2 + 3.x - 4 = 0   and   2.x2 + 3.x + 4 = 0

             2  ENTER^
             3  ENTER^
           -4  XEQ "P2"  >>>>   -2.350781059
                    X<>Y                   0.850781060      Flag  F00 is clear, therefore the 2 solutions are  -2.350781059  and  0.850781060

            2  ENTER^
            3  ENTER^
            4     R/S         >>>>   -0.75
                 X<>Y                    1.198957881      Flag  F00 is set, therefore the 2 solutions are  -0.75 + 1.198957881.i  and   -0.75 - 1.198957881.i
 
 

Cubic Equation
 

"P3" finds the 3 roots of   a.x3+b.x2+c.x+d  by the Cardan's ( or Tartaglia's ) formulae  ( with a # 0 )
 
 
 
      STACK        INPUTS      OUTPUTS
           T          a # 0             /
           Z             b             z
           Y             c             y
           X             d             x

-If  CF 01  the 3 solutions are  x ; y ; z
-If  SF 01   ------------------  x ; y+i.z ; y-i.z

Example:      Solve  2.x3-5.x2-7.x+1 = 0  and   2.x3-5.x2+7.x+1 = 0

        2   ENTER^
       -5  ENTER^
       -7  ENTER^
        1   XEQ "P3" >>>>   3.467727065   RDN   0.131206041   RDN   -1.098933107   which are the 3 real solutions since flag F01 is clear.

        2  ENTER^
       -5  ENTER^
        7   ENTER^
        1     R/S         >>>>   -0.130131632   RDN   1.315065816   RDN   1.453569820

-Flag F01 is set , therefore the 3 solutions are  -0.130131633  ;  1.315065816 - 1.453569820.i  ;  1.315065816 + 1.453569820.i
 
 

Quartic Equation
 

"P4" solves the equation   x4+a.x3+b.x2+c.x+d = 0   ( if the leading coefficient is not 1 , divide all the equation by this coefficient )
 
 
 
      STACK        INPUTS      OUTPUTS
           T             a             t
           Z             b             z
           Y             c             y
           X             d             x

-If  CF 01  &  CF 02  the 4 solutions are    x ; y ; z ; t
-If  SF 01  &  CF 02   ------------------    x+i.y ; x-i.y ; z ; t
-If  SF 01  &  SF 02   ------------------    x+i.y ; x-i.y ; z+i.t ; z-i.t

Example1:    Solve  x4-2.x3-35.x2+36.x+180 = 0

       -2   ENTER^
      -35  ENTER^
       36   ENTER^
      180  XEQ "P4"  >>>>   6  RDN   3   RDN   -2   RDN   -5    Flags F01 & F02 are clear , the 4 solutions are   6 ; 3 ; -2 ; -5

Example2:    Solve  x4-5.x3+11.x2-189.x+522 = 0

        -5   ENTER^
        11   ENTER^
     -189   ENTER^
       522      R/S     >>>>   -2  RDN    5   RDN   3   RDN   6  Flag F01 is set & F02 is clear , the 4 solutions are  -2+5.i ; -2-5.i ; 3 ; 6

Example3:    Solve  x4-8.x3+26.x2-168.x+1305 = 0

       -8   ENTER^
       26   ENTER^
     -168  ENTER^
     1305    R/S      >>>>   -2  RDN    5   RDN   6   RDN   3    Flags F01 & F02 are set , the 4 solutions are  -2+5.i ; -2-5.i ; 6+3.i ; 6-3.i
 

Carlson Elliptic Integral 1st kind: RF(x,y+iz,y-iz)
 
 
 
 
      STACK       INPUTS    OUTPUTS
          Z            z           /
          Y            y           /
          X            x  RF(x,y+iz,y-iz)

 
Example:

     4  ENTER^
     2  ENTER^
     1  XEQ "RFZ"  >>>>  RF ( 1 , 2+4i , 2-4i ) = 0.631488738
 

Carlson Elliptic Integral 3rd kind: RJ(x,y+iz,y-iz,p)
 

 RJZ  computes  RJ ( x , y+i.z , y-i.z , p )  with  p > 0
 
 
        STACK        INPUTS     OUTPUTS
            T           p>0            /
            Z             z            /
            Y             y            /
            X             x   RJ(x,y+iz,y-iz,p)

 
Example:

    4  ENTER^
    3  ENTER^
    2  ENTER^
    1  XEQ "RJZ"  >>>>  RJ ( 1 , 2+3i , 2-3i , 4 ) = 0.205644141
 
 
 

References:
 

 [1]   Stamatia Mavridès - "L'Univers relativiste" - Masson  ISBN 2-225-36080-7  ( in French )
 [2]   Jean Heidmann - "Introduction à la cosmologie" - PUF  ( in French )
 [3]   Landau & Lifshitz - "Classical Theory of Fields" - Pergamon Press
 [4]   Feige - "Elliptic Integrals for Cosmological Constant Cosmologies" - Astron. Nachr 313 (1992) 3, 139-163
 [5]   Kantowski, Kao, Thomas - "Distance-Redshift Relations in Inhomogeneous Friedmann-Lemaitre-Robertson-Walker Cosmology" -
        The Astrophysical Journal, 545: 549-560
 [6]   Michèle Caputo - "Some Space Gravity Formulas and the Dimensions and the Mass of the Earth"
 [7]   http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
 [8]   Fairhead & Bretagnon - "An Analytica Formula for the Time Transformation TB-TT" Astronomy & Astrophysics 229, 240-247 ( 1990 )
 [9]   http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
[10]  Henri Arzelies - "Relativité Généralisée, Gravitation" - Gauthier-Villars ( in French )
[11]  Vincenty's formula
[12]  Paul D. Thomas - "Spheroidal Geodesics, Reference Systems & Local Geometry" - US Naval Oceanographic Office.
[13]  Marcin Ligas  - Cartesian to geodetic coordinates conversion on a triaxial ellipsoid - J Geod (2012) 86:249–256
[14]  If you want to visualize the geodesics on a triaxial ellipsoid, go to this excellent page
[15]  Nachum Dershowitz & Edward M. Reingold - "Calendrical Calculations" - Cambridge University Press - ISBN 978-0-521-70238-6
[16]  Nachum Dershowitz & Edward M. Reingold - "Indian Calendrical Calculations"