Orbit Determination for the HP48G/GX&HP49/50G

-This zip file contains:

1-A binary file ORB48G ( #46744d / 16185 bytes ) for the HP48G/GX
2-A binary file ORB50G ( #13146d / 16328.5 bytes ) for the HP50G

3-A text file ORB50 with the same programs as ORB50G

4-A 3x6 matrix for the example in GAUSSHG in the txt file DATA1
5-A 6x6 matrix for the example in HERGET in the txt file DATA2

6-This html page

-ORB48G & ORB50G are directories.
 

Notes:

-In all the examples, I've used an HP48GX.
-Though the programs are almost identical, the results given by an HP50
  may be slightly different from the results given by an HP48.
-Execution times are of course smaller with an HP50.
 

Warning:

-Flags 1 , 2 , 4 are used in these programs.
-There are indicators for these flags on the HP48 but not on the HP50
-So, check the state of these flags carefully if you have an HP50 !

 1  CF  =  Standard Ecliptic  J2000.0 = 2000/01/01 12h TT
 1  SF  =   Ecliptic of the Date

 4  CF  = Gregorian Calendar
 4  SF  =  Julian Calendar

-Flag 2 is also used in the p-iteration method  ( cf  RR->V )

  2  CF =  Angle ( Comet Position1 - Sun - Comet Position2 )  < 180°
  2  SF =  Angle ( Comet Position1 - Sun - Comet Position2 )  > 180°
 
 
 

Functions

Description

GAUSSHG    SUN      DT->T      J2000      J->DAT      T      PRE0      SUN      SP       SCOS  TADXYZ  TLBXYZ  HERGET  IMPROVE  EPS  RR->V  R  V  RV->EL  RV->R  T  q  E  I  w  W  n  p  a  IwW  PHEM  VP  KPL  S->R  R->S  EE  AD->LB  OBL  P->R  R->P  KG  OUT

Gauss-Herrick-Gibbs Preliminary Orbit Determination    Directory = Sun      Date & Time -> T      Date -> Number of days since 2000/01/01 12h TT      Number of days -> Date      Time in millenia since J2000.0      Correction for precession      Position of the Sun      Sum of Poisson series      Sum of trigonometric series  DATA = Time-Rightascension-Declination-Rectangular Coordinates  DATA = Time-Longitude-Latitude-Rectangular Coordinates  Orbit Determination by the method of Herget ( >3 observations )  Improving an Orbit  Estimation of the error  2 positions -> Velocity1 ( p-iteration method )  Vector-Position  Vector-Velocity  Position+Velocity -> Orbital Element  Position1+Velocity1 -> Position2  Instant of passage in Perihelion  Perihelion Distance  Eccentricity  Inclination  Argument of the Perihelion  Longitude of the Ascending Node  Mean Motion  Parameter  Semi-major Axis  Reduction of Ecliptic Elements from one Equinox to another one  Parabolic-Hyperbolic-Elliptic Motion  M , e#1  --->  True Anomaly v & r/q  Kepler's Equation:   M , e#1  --->  Eccentric anomaly E  Spherical->Rectangular Conversion  Rectangular->Spherical Conversion  Equatorial<->Ecliptic Conversions  Takes TADXYZ and creates TLBXYZ  Obliquity of the Ecliptic  Polar->Rectangular Conversion  Rectangular->Polar Conversion  Gauss Gravitational Constant = 0.01720209895  Places the orbital elements T q E I w W  in level 6-5-4-3-2-1

 

Gauss-Herrick-Gibbs Preliminary Orbit Determination
 

-Given 3 observations of an asteroid, comet... "GAUSSHG" returns the orbital elements of its orbit.

-Gauss' method computes approximate values of  T , q , e , i , omega , OMEGA  from 3 observations of a comet.
-The method fails if the inclination i = 0 or is very small.

-With the Herrich-Gibbs refinement, the method becomes 4th order in time !

-At an instant  t₁  the right ascension of the comet =  RA₁  its declination = d₁ and the rectangular ecliptic coordinates of the Sun are  X₁  Y₁  Z₁
-At an instant  t₂  the right ascension of the comet =  RA₂  its declination = d₂ and the rectangular ecliptic coordinates of the Sun are  X₂  Y₂  Z₂
-At an instant  t₃  the right ascension of the comet =  RA₃  its declination = d₃ and the rectangular ecliptic coordinates of the Sun are  X₃  Y₃  Z₃

-We must have   t₁ <  t₂ <  t₃

 >>> The results will be more accurate if  t₃ - t₂ is equal ( or nearly equal ) to   t₂ - t₁  ( this value is often of the order of a few days )

-Let

           l_i  =  cos RA_i cos d_i
           m_i = sin obl sin d_i + cos obl cos d_i sin RA_i    where   obl is the obliquity of the ecliptic,   i = 1 , 2 , 3
           n_i = cos obl sin d_i - sin obl cos d_i sin RA_i

-Let D_i = distance Earth-Comet at the instant t_i , u_i ( l_i , m_i , n_i ) , -R_i = position vector of the Earth , r_i = position vector of the Comet

    r_i = D_i u_i  - R_i        (  || u_i || = 1  )

-If we neglect the planetary perturbations, the 3 position vectors r_i  lie in the same plane, so there are coefficients c₁ , c₃ such that

   r₂ = c₁ r₁ + c₃ r₃   and likewise for the velocity at the second instant:        V₂ = -d₁ r₁ + d₂ r₂ + d₃ r₃

-One can prove that, approximately:     c_i = A_i + B_i / r₂³ ( 1 + B₂ / r₂³ )   ( i = 1 , 3 )  and   d_i = G_i + H_i / r_i³  ( i = 1 , 2 , 3 )

     where    A₁ = (t₃ - t₂)/(t₃ - t₁)  ,   B₁ = k²A₁ [ (t₃ - t₁)² - (t₃ - t₂)² ]/6      with   k =  0.01720209895 = Gaussian gravitational constant
                  A₃ = (t₂ - t₁)/(t₃ - t₁)  ,   B₃ = k²A₃ [ (t₃ - t₁)² - (t₂ - t₁)² ]/6

      and      G₁ = (t₃ - t₂)/[(t₃ - t₁)(t₂ - t₁)]       G₃ = (t₂ - t₁)/[(t₃ - t₁)(t₃ - t₂)]        G₂ = G₁ - G₃
                 H₁ = k²(t₃ - t₂)/12                        H₃ = k²(t₂ - t₁)/12                         H₂ = H₁ - H₃

-Let   A( A₁ X₁ - X₂ + A₃ X₃ , A₁ Y₁ - Y₂ + A₃ Y₃ ,  A₁ Z₁ - Z₂ + A₃ Z₃ )
 and   B( B₁ X₁ + B₃ X₃ , B₁ Y₁ + B₃ Y₃ ,  B₁ Z₁ + B₃ Z₃ )

         P₁ = A.( u₂ x u₃ )/D          Q₁ = B.( u₂ x u₃ )/D
         P₂ = A.( u₁ x u₃ )/D          Q₂ = B.( u₁ x u₃ )/D          where      D = u₁.( u₂ x u₃ )           ( . = dot product , x = cross product )
         P₃ = A.( u₁ x u₂ )/D          Q₃ = B.( u₁ x u₂ )/D

-The distance Earth-Comet at the central instant is found by solving the equation

        D₂ = P₂ + Q₂ / r₂³ ( 1 + B₂ / r₂³ )  with    r₂ =  D₂² - 2 D₂ ( l₂ X₂ + m₂ Y₂ + n₂ Z₂ ) + X₂² + Y₂² + Z₂²

-Then,  D₁ = ( P₁ + Q₁ / r₂³ ( 1 + B₂ / r₂³ ) )/( A₁ + B₁ / r₂³ ( 1 + B₂ / r₂³ ) )
  and    D₃ = ( P₃ + Q₃ / r₂³ ( 1 + B₂ / r₂³ ) )/( A₃ + B₃ / r₂³ ( 1 + B₂ / r₂³ ) )

  and the rectangular ecliptic coordinates of the comet at the 3 instants are:

           x_i = l_i D_i - X_i
           y_i = m_i D_i - Y_i        i = 1 , 2 , 3         so we know the 3 vectors r_i's   and the velocity V₂ = -d₁ r₁ + d₂r₂ + d₃ r₃   may be computed
           z_i = n_i D_i - Z_i

-So the position and velocity vectors at the second instant are computed and we can call RV->EL
 

>>>>  First store the required data in 'TADXYZ' and then, execute  GAUSSHG
 
 

STACK	INPUTS	OUTPUTS
Level 7	/	D2
Level 6	/	T
Level 5	/	q
Level 4	/	E
Level 3	/	I
Level 2	/	w
Level 1	D	W

   Where  D  is an estimtion of the distance Sun-Comet at the 2nd instant and D2 the computed distance Earth-Comet at this instant.

Example:        Comet C/2014 AA52 Catalina, on 2015/02/01 0h TT , 2015/02/10 0h  TT , 2015/02/20 0h , astrometric coordinates ( /J2000 ):

                   2015/02/01 0h                   2015/02/10 0h                  2015/02/20 0h

      t                 5509.5                              5518.5                              5528.5                               days
    RA        1h07m43s058                     0h58m40s151                   0h53m53s415
   Decl       -57°17'23"42                      -52°05'21"91                    -46°54'15"67
     X          0.653892160                      0.763553245                   0.863088915                  Astronomical Units      X  Y  Z  are given by the SOLEX
     Y        -0.736974521                     -0.624900515                 -0.482202751                   Astronomical Units      software ( cf reference [1] )
     Z          0.000019390                      0.000019018                   0.000014378                   Astronomical Units
 

-These data are the first 3 rows of the matrix in 'TADXYZ': each row = t  RA Decl  X Y Z )

>>> So, the columns above = the rows in 'TADXYZ'

-Store the matrix DATA1 in the variable 'TADXYZ'

[[ 5509.5 1.0743058 -57.172342 .65389216 -.736974521 .000019390 ]
 [ 5518.5 .5840151 -52.052191 .763553245 -.624900515 .000019018 ]
 [ 5528.5 .5353415 -46.541567 .863088915 -.482202751 .000014378 ]]
 

-If you choose  D = 3 AU

  1 CF   3   GAUSS   >>>>   D2 = 2.403691                                            ---Execution time = 16s---     with an HP48GX

                                               T = 5536.18831 = 2015/02/27.6883 TT
                                               q  =   2.002315
                                               e  =   0.999457
                                               i   = 105°2133
                                      omega   = -67°7290
                                   OMEGA = -29°5133

  and we get in the variables 'n'  'P'  'a'

                                           n      =  0.000004407 deg/d
                                           P     =  4.003542
                                           a      =  3684.26703

Notes:

-During the calculation, the HP-48 solves an equation with the nuilt-in function  ROOT
-If a negative solution is found for D, the program stops and displays:  Error: Bad Guess(es)
-Place another value in level 1 and press GAUSSHG again

-In some cases, the problem may have 2 or even 3 possible solutions
-Other observations are required to choose the correct solution.

-Remember that R & V are calculated at the second instant of the DATA in 'TADXYZ'

-The algorithm is very sensitive to roudoff-errors
 

Time in millenia since J2000.0
 

 DT->T  takes the date in level 2 and the time in level 1
 and calculates the time in millenia since J2000.0 i-e 2000/01/01  12h TT
 
 

STACK	INPUTS	OUTPUTS
Level 2	YYYY.MNDD	/
Level 1	HH.MNSS	/

Example:       On  2015/05/31  at  12h41m24s77  TT

              4           CF
      2015.0531    ENTER
      12.412477     DT->T    stores   1.54114408184 E-2  in  'T'
 

Warning:

-If flag 4 is clear, the date are Gregorian dates
-If flag 4 is set, the dates are Julian dates.
 

Date -> Number of days since 2000/01/01 12h TT
 
 
 

STACK	INPUT	OUTPUT
Level 1	YYYY.MNDD	days

Example:       On  2015/05/31  at  18 h

  4  CF     2015.053175   J2000   returns   5629.25
 

Warning:

-If flag 4 is clear, the date are Gregorian dates
-If flag 4 is set, the dates are Julian dates.
 

Number of days -> Date
 

  J->DAT  performs the reverse operation
 
 

STACK	INPUT	OUTPUT
Level 1	days	YYYY.MNDD

Example:     5629.25   J->DAT   >>>>   2015.053175
 

Warning:

-If flag 4 is clear, the date are Gregorian dates
-If flag 4 is set, the dates are Julian dates.
 

Correction for precession
 

 PRE0  takes the ecliptic longitude & latitude at J2000.0 ( i-e on 2000/01/01 12h TT )
 and returns the ecliptic longitude & latitude of the date specifyed by the value in 'T' ( in millenia since J2000.0 )
 
 
 

STACK	INPUTS	OUTPUTS
Level 2	L2000	L
Level 1	B2000	B

    Where all the angles are in decimal degrees

Example:     The eclptic coordinates at the standard epoch J2000.0  are  L₂₀₀₀  = 49°7891   B₂₀₀₀ = 24°1279

-Calculate L & B for T = 1000 julian years since J2000 ( i-e  3000/01/08 12 TT )

 1  'T'  STO

 49.7891    ENTER
 24.1279     PRE0             >>>>      L = 63°757208
                                                        B = 24°236957
 

Position of the Sun
 

-This program is not very useful if you employ a very accurate software like SOLEX ( cf reference [1] )
-Otherwise, SUN will give you the ecliptic rectangular coordinates of the Sun with an error smaller than E-6 AU for the time span [1900-2100]
  and still a good precision several millenia around J2000.

   1 CF  =  Coordinate with respect to the standard ecliptic J2000
   1  SF =  ------------------------------- ecliptic of the date
 
 

STACK	INPUTS	OUTPUTS
Level 3	/	X
Level 2	/	Y
Level 1	/	Z

   Store T in the variable 'T' before executing SUN ,  X   Y   Z  are expressed in Astronomical Units  ( 1 AU = 149597870.7 km )

Example:    On 2100/04/04  21h41m37s  TT

   2100.0404    ENTER
     21.4137      DT->T             ( T = 0.100255725943 is now stored in 'T' )

  •  Standard ecliptic J2000

     1 CF     SUN    gives          X = 0.9710287  AU                   ---Execution time = 8.0s---   with an HP48GX
                                               Y = 0.2386091  AU
                                               Z = -0.0000776  AU

  •  Ecliptic of the Date

     1 SF     SUN    gives          X = 0.9649056  AU                   ---Execution time = 8.4s---   with an HP48GX
                                               Y = 0.2622752  AU
                                               Z = -0.0000026  AU

Warning:

-If flag 4 is clear, the date are Gregorian dates
-If flag 4 is set, the dates are Julian dates.
 

DATA = Time-Right-ascension-Declination-Rectangular Coordinates
 

'TADXYZ'  must contain a nx6 matrix ( n = number of observations )

 [ [  t1   RA1  Decl1  X1  Y1  Z1  ]
   [  t2   RA2  Decl2  X2 Y2  Z2  ]
   [  t3   RA3  Decl3  X3  Y3  Z3  ]                at least 3 observations
   ..................................................

   [  tn   RAn  Decln  Xn  Yn  Zn  ] ]

Where  t  must be expressed in days since  J2000.0 = 2000/01/01  12h TT
           RA --------------------  hh.mnss
         Decl --------------------   sexagesimal degrees  ° ' ""
         X Y Z = rectangular ecliptic coordinates of the Sun expressed in Astronomical Units ( 1AU = 149597870.7 km )
 

DATA = Time-Longitude-Latitude-Rectangular Coordinates
 

The DATA in 'TADXYZ' are automatically transformed in 'TLBXYZ' by AD->LB
 
 

Orbit Determination by the method of Herget ( 3 or more than 3 observations )
 

-Herget's method regards the 1st and last observations as exact.

-Assuming there are n observations, you have to give an estimate of the distances Earth-Asteroid at the first and last instants  D₁ & D_n
 and an estimate of the parameter p.

-So, the HP-48 can compute the position-vectors r_i  for  i = 2 , 3 , .... , n-1  thanks to "RR-V" & "RV-R"

-For i = 2 , .... , n-1  the dot-products

   P_i = ( r_i + R_i ).A_i               where      R_i = rectangular ecliptic coordinates of the Sun                                                  L_i = eclptic longitude
   Q_i = ( r_i + R_i ).D_i                          A_i = ( -Sin L_i , Cos L_i , 0 ) D_i = ( -Sin b_i Cos L_i , -Sin b_i Sin L_i , Cos b_i )          b_i = ecliptic latitude

 ( Herget uses the right-ascensions & declinations instead of the longitudes & latitudes )

-If there were no errors, we'd have  P_i = Q_i = 0
-To find the correct D₁ & D_n , we solve non-linear systems of  ( 2n-4 ) equations in 2 unknowns by Newton's method & least-squares approximation.

    (¶P_i/¶D₁)DD₁ + (¶P_i/¶D_n)DD_n = - P_i
    (¶Q_i/¶D₁)DD₁ + (¶Q_i/¶D_n)DD_n = - Q_i

-The partial derivatives are approximated by formulas of the type: ¶f/¶x ~ [ f(x+h) - f(x) ] / h

Warning:

    Clear F02 if the angle  ( r₁ ; r_n ) < 180°
     Set   F02 if the angle  ( r₁ ; r_n ) > 180°
 
 

STACK	INPUTS	OUTPUTS1	OUTPUTS2	OUTPUTS3 ...	OUTPUT
Level 4	/	/	D'1	D"1	orbital elements
Level 3	h	/	D'n	D"n	orbital elements
Level 2	D1	D'1	D"1	D"'1  ...	orbital elements
Level 1	Dn	D'n	D"n	D"'n   ...	orbital elements

   Where  D₁ & D_n  are 1st guesses of the distance Earth-Comet at the first & last instant ( D is denoted RHO on the stack )
                    h = a small value to approximate the partial derivatives  ( usually  0.01 AU  works well )

    HERGET also uses the value in the variable 'P' as a 1st guess of the parameter

Example:   Let's take again Comet C/2014 AA52 Catalina, astrometric coordinates ( /J2000 ) with the 6 observations below,

  where the right-ascensions have been rounded to 0.1 second and the declinations to 1 arcsecond
 

      t                 5509.5                    5518.5                      5528.5                    5537.5                 5546.5                     5556.5
    RA          1h07m43s1              0h58m40s2                0h53m53s4            0h52m18s7          0h52m13s9             0h53m10s1
   Decl         -57°17'23"               -52°05'22"                -46°54'16"              -42°45'51"          -39°04'47"               -35°27'29"
     X          0.653892160           0.763553245            0.863088915          0.930110731       0.974274312           0.995480570
     Y        -0.736974521           -0.624900515          -0.482202751       -0.341009813      -0.191511107         -0.020002821
     Z          0.000019390            0.000019018            0.000014378         0.000005490        0.000004634         -0.000001613
 

1st Step      Store the matrix in DATA2 in the variable  'TADXYZ'
 

[[ 5509.5 1.07431 -57.1723 .65389216 -.736974521 .00001939 ]
 [ 5518.5 .58402 -52.0522 .763553245 -.624900515 .000019018 ]
 [ 5528.5 .53534 -46.5416 .863088915 -.482202751 .000014378 ]
 [ 5537.5 .52187 -42.4551 .930110731 -.341009813 .00000549 ]
 [ 5546.5 .52139 -39.0447 .974274312 -.191511107 .000004634 ]
 [ 5556.5 .53101 -35.2729 .99548057 -.020002821 -.000001613 ]]
 

2nd Step      Before using "HERGET" , let's execute "GAUSSHG" with the first 3 observations

-With  D = 3 in Level 1 ,   GAUSSHG  yields:              D2 = 2.406895  AU

                                               T = 5535.913678 = 2015/02/27.4137 TT
                                               q  =   2.005546
                                               e  =   1.007843
                                               i   = 105°1904
                                      omega   = -67°8108
                                   OMEGA = -29°4423

  and we get in the variable  'P'   P = 4.026821 AU

3rd Step     Now, we are ready to use 'HERGET'

-A good choice for h is often  0.01  AU
-The results given by "GAUSSHG"   suggest    D₁ =  D_n = 2.407  thus,

                           2  CF

         0.01           ENTER
        2.407          ENTER   HERGET  >>>>    RHO1' = 2.324644316               ---Execution time = 2m45s---       ( with an HP48GX )
                                                                        RHON' = 2.724885559

-To get the next approximation(s), press  left shift CONT without clearing the stack ,  it yields

                        CONT  >>>>   RHO1" = 2.315488110
                                                 RHON" = 2.715577781

                        CONT  >>>>   RHO1"' = 2.315402260
                                                 RHON"' = 2.715477404

                        CONT  >>>>   RHO1"" = 2.315401995
                                                 RHON"" = 2.715477097

-Here, the new position- & velocity- vectors have not yet been calculated.
 

4th Step       If you want to stop the iterations and get the orbital elements, clear the stack and press  left shift  CONT

             CLEAR   CONT   >>>>          T = 5536.16433 = 2015/02/27.6643

 and               q      =   2.002984
                      e      =   1.000995
                      I      = 105°20839
                 omega   = -67°7296  = 292°2704
               OMEGA = -29°4981  = 330°5109

Notes:

-The last D-values are saved in level 8 and level 7
-The method would also converge with D₁ =  D_n = 2 , but better guesses will reduce the number of iterations.
-Remember that R & V are calculated at the first instant of the DATA in 'TADXYZ'

-Of course, the execution time increases with the number of observations
 

Comparison with the exact elements:

-The exact orbital elements - given for instance by the imcce - are

              T =  2015 feb 27,64787 TT ±0,00000 TT = 5536.14787  TT
              q =  2,0025966  ±0,0000008
               e =  1,0004430  ±0,0000019
               i =  105,2112331°  ±0,0000262°
      omega =  292,2632213°  ±0,0000207°
  OMEGA = 330,4930204°  ±0,0000259°

-The JPL gives almost the same results:

 Element                  Value                                    Uncertainty (1-sigma)    Units

     T           2015-Feb-27.61499079                         0.00040407
     q             2.002902188984628                             3.7736e-06              AU
     e             1.000563179802131                             6.0681e-06
      i              105.207183638451                               3.7436e-05             deg
   peri           292.2449325905247                             0.00018586             deg
  node          330.4895894198529                              2.3831e-05             deg
 

Improving an Orbit
 

-After finding a preliminary orbit determination, you might want to improve your results.

-IMPROVE  uses the values in 'R' & 'V' at a given time t
 to evaluate the variation in the positions that results if a coordinate of R , say x , is replaced by x+h
 and similarly if a coordinate of V , say  V_x  , is replaced by   V_x + h'  ( this program uses h' = h/100 )

-The partial derivatives are approximated by formulas of the type  df/dx ~ [ f(x+h) - f(x) ] / h

-This leads to an overdetermined linear system which is solved by a least-square method, via the built-in function  LSQ

-Gradually,  R & V  are improved until you judge that the changes become negligible.
-Finally, RV->EL is called to give the orbital elements.
 
 

STACK	INPUT	OUTPUT	FINAL OUTPUT
Level 2	t	/	orbital elements
Level 1	h	[dx dy dz dvx dvy dvz ]	orbital elements

   Where  t  is the number of days since J2000 for which 'R' & 'V'  have been calculated
               h is a small increment to approximate the partial derivatives ( usually h = 0.01 AU )
      and   [dx dy dz dvx dvy dvz ]  is the last correction to the vectors R & V

Example:   If we start with the results given by GAUSSHG at the 2nd step of HERGET

   5518.5     ENTER
     0.01       IMPROVE    >>>>   [ -1.511 E-3  4.297 E-4  1.990 E-3  -2.611 E-5  -4.082 E-6  -5.187 E-6 ]     ---Execution time = 131s---  ( HP48GX )

-Press left shiftt  CONT  without clearing the stack  to get the next improvement(s)

            CONT   >>>>   [ -2.970 E-5  1.218 E-5  5.246 E-5  -6.838 E-7  -1.787 E-7  -8.764 E-8 ]

            CONT   >>>>   [ -5.093 E-7  2.251 E-7  9.925 E-7  -1.329 E-8  -3.266 E-9  -1.025 E-9 ]

-If you think the next improvement would be negligible, clear the stack and CONT

    CLEAR   CONT   >>>>             T = 5536.17285 = 2015/02/27.6729

 and               q      =   2.002883
                      e      =   1.000748
                      I      = 105°20911
                 omega   = -67°7271  = 292°2729
               OMEGA = -29°5004  = 330°4996

Note:

-Of course, the execution time increases with the number of observations
 

Estimation of the error
 

-EPS calculates the square-root of the sum of the squares of the angular separations
 between the observed positions and the calculated positions ( in decimal degrees ).

-We assume that 'R' & 'V' contain the position- & velocity- vectors at the instant ( in days since J2000.0 ) specified in level 1
  and that 'TLBXYZ' contains the required data.
 
 

STACK	INPUT	OUTPUT
Level 1	t	eps

   Where  t  is the number of days since J2000 for which 'R' & 'V'  have been calculated

Examples:

1-With the result found by GAUSSHG at the 2nd step of HERGET,

      5518.5   EPS  gives   2.174 E-3 deg = 7"827

2-With the results finally given by HERGET,

      5509.5   EPS  returns  2.581 E-4 deg = 0"929

3-With  IMPROVE  in the paragraph above,

      5518.5   EPS  yields  2.538 E-4 deg = 0"914   only slightly better than Herget's results

4-With the 1st orbital elements at the top of this page with GAUSSHG - which were found with much more accurate data -
    but compared to the 6 observations in DATA2

     5518.5   EPS  >>>  5.699 E-4 deg = 2"052
 

Position1+Velocity1 -> Position2
 

-We could use  RV->EL and  PHEM  but in order to avoid a possible loss of accuracy if the eccentricity is very close to 1,
 "RV-R" employs the universal variable x  to perform the calculations ( cf reference [2] )

-We solve the equation:   f(x) = x³ S(z) + ( r.V / k ) x² C(z) + r x ( 1 - z S(z) ) -  t k = 0   by the method of Newton

   here,  z = x² / a    where   a = semi-major axis              1 / a = 2 / r - V² / k²

       t = ( t₂ - t₁ )  and  k = 0.01720209895 = Gauss gravitational constant

  C(z) = ( 1 - Cos z^1/2 ) / z = ( 1 - Cosh (-z)^1/2 ) / z                  S(z) =  ( z^1/2 - Sin (z^1/2) ) / ( z^3/2) = ( Sinh (-z)^1/2 - (-z)^1/2 ) / (-z)^3/2

-After founding  x  ,  we have

      f = 1 - x² C(z) / r
      g = t - x³ S(z) / k

  and   r' = f r + g V
 
 

STACK	INPUT	OUTPUT
Level 1	t2-t1	R2

Example:        R₁ = [ 0.16 , 1.38 , 0.24 ]    V₁ = [ 0.015 , 0.01 , 0.001 ]    Find the position 100 days later

-Store these 2 vectors in the variables 'R' & 'V' respectively

     100  RV->R  returns   R₂ = [ 1.50929963708 , 1.91954203011 , 0.265117222709 ]   in level 1
 

2 positions -> Velocity
 

-We solve Lambert problem: given 2 position-vectors  r₁ & r₂  at different instants  t₁ & t₂
  "RR-V"  stores the velocity vector  V₁ at the first instant in the variable 'V'

-It uses an iteration method called the method of Herrick & Liu  or  p-iteration method ( cf reference [3] )

>>> Flag F02 must be cleared if the angle  ( r₁ ; r₂ ) < 180°
>>> Flag F02 must be set if the angle  ( r₁ ; r₂ ) > 180°

-Let    C = Cos ( v₂ - v₁ ) = ( r₁.r₂ ) / ( r₁. r₂ )          where    v  is the true anomaly                  Sin ( v₂ - v₁ ) = +/- ( 1 - C² )^1/2

-We have   f = ( r₂ / p ) ( C - 1 ) + 1   &  g = ( r₁. r₂ ) / ( k p^1/2 )  Sin ( v₂ - v₁ )

  from which   V₁ =  ( r₂ - f r₁ ) / g   whence  an estimate of r₂ ( say r* )   whence   C* = ( r₁.r* ) / ( r₁. r* )

-The solution is found if  C - C* = 0
 
 

STACK	INPUTS	OUTPUTS
Level 4	t1	/
Level 3	R1	/
Level 2	t2	/
Level 1	R2	/

   When the program stops,  the velocity vector V₁ at the first instant in the variable 'V'

Example:    t₁ = 0   ,   R₁ = [ 0.16 , 1.38 , 0.24 ]  ,   t₂ = 100 days  ,  R₂ = [ 1.50929963708 , 1.91954203011 , 0.265117222709 ]

-If you choose  2  CF  and  p = 1  ,

-Store  1  in the variable 'P'

                                            2                                               CF
                                            0                                           ENTER
                             [ 0.16 , 1.38 , 0.24 ]                             ENTER
                                          100                                          ENTER
  [ 1.50929963708 , 1.91954203011 , 0.26511722709 ]   RR-V  >>>>         /                            ---Execution time = 12s---  with an HP48GX

 and the velocity-vector in 'V' = [ 0.015 , 0.01 , 0.001 ]     ( with small differences in the last decimals )
 

Notes:

 R₁ is saved in the variable 'R'
 The parameter p is in the variable 'p'  here,  p =1.27633801262

-If the 1st guess for P is too bad, it will lead to an error while solving the equation,
-In this case, store another - usually smaller - value in 'P' and start again.
 

Position+Velocity -> Orbital Element
 

-If we know the position and velocity vectors R & V  at an instant t , we can deduce the 6 orbital elements  T  q  e  i  omega  OMEGA
 

-We calculate the cross-product  h = r₂ x V₂ and then:

-The parameter  p = h²/k²
-The inclination  i = Acos h_z/h = Atan2 (h_x²+h_y²)^1/2/h_z

-We have also  V² / k² - 2 / r = 1 / a   whence the eccentricity e

     r = p / ( 1 + e Cos v )

-The longitude of the ascending node  OMEGA = Atan2 h_x/(-h_y)

  ( v + omega ) may be computed by the formula:

    v + omega = sign(z) Acos [ ( x Cos OMEGA + y Sin OMEGA ) / r ]

-And the time of passage in perihelion is computed via Kepler's equation.
 
 

STACK	INPUT	OUTPUT
Level 1	t	/

  The orbital elements are in the variables  T q E I w W n p a

Example:      With the vectors R & V above:  R = [ 0.16 , 1.38 , 0.24 ]  , V = [ 0.015 , 0.01 , 0.001 ]  and assuming  t = 0

   0   RV-EL  >>>>  /

 and we find in the corresponding variables ( execute OUT to place them on the satck )

                      T     =  -76.072511 days
                      q      =   0.720414  AU
                      e      =   0.771672
                       i      = 169°36068
                 omega   =  15°72463
               OMEGA = -163°48430

 moreover:
                      n      =   0.175861 deg/d
                      p     =   1.276338  AU
                      a      =  3.155171  AU
 

Reduction of Ecliptic Elements from one Equinox to another one
 

-The 3 following orbital elements of a planet or a comet change when they are referred to different equinoxes

    i =  inclination
   w = argument of perihelion
   W = longitude of ascending node
 
 

STACK	INPUTS	OUTPUTS
Level 5	i1( deg )
Level 4	w1 ( deg )
Level 3	W1  ( deg )	i2 ( deg )
Level 2	YYYY.MNDDdd1	w2 ( deg )
Level 1	YYYY.MNDDdd2	W2  ( deg )

    where   i₁ w₁W₁ are reffered to the 1st date ( level2 )   and    i₂ w₂W₂ are reffered to the 2nd date ( level1 ) , all expressed in decimal degrees

Example:

-The elements of a comet are    i₁ = 47°1234  ,  w₁ = 151°5432  , W₁ = 34°5678   on  1984/12/28
-Calculate these elements on  2345/07/24

      47.1234   ENTER
     151.5432  ENTER
      34.5678   ENTER
  1984.1231   ENTER
  2345.0724     iwW        returns

     i₂  =  47°1590
    w₂ = 151°5850
   W₂ =  39°5796
 

Parabolic-Hyperbolic-Elliptic Motion
 

-Assuming the orbital elements  T q E I w W  are stored in the corresponding variables,
  PHEM  takes t ( in days since J2000.0 ) in level 1
  and returns [ X Y Z ] the rectangular ecliptic coordinates in level 1
 
 

STACK	INPUT	OUTPUT
Level 1	t	/

  Assuming the orbital elements are in the variables  T q E I w W

Example:      With the elements found in the example RV->EL

     100  PHEM  gives  [ 1.50929963708 , 1.91954203009 , 0.265117222712 ]

-Almost the same result we found with  RV->R
 

Kepler's Equation
 

KPL takes the mean anomaly M in level 2 and the eccentricity e # 1 in level 1
It returns the eccentic anomaly E
 
 

STACK	INPUTS	OUTPUTS
Level 2	M ( deg )	/
Level 1	e # 1	E ( deg )

Examples:

    0.4    ENTER
    1.2     KPL        gives    E = 0.987853766458

     7     ENTER
    0.9     KPL        yields   E = 40°4669345175
 

Spherical -> Rectangular Conversion
 

-This program employs the usual formulae ( for the astronomers ):

       X = R cos B cos L
       Y = R cos B sin L
       Z = R sin B
 
 

STACK	INPUTS	OUTPUTS
Level 3	R	X
Level 2	L	Y
Level 1	B	Z

Example:     R = 7   L = 124°   B = 41°

                        7   ENTER
                      124 ENTER
                       41  S->R      >>>>   X = -2.95419769010
                                                       Y =  4.37977818860
                                                       Z =  4.59241320294

Notes:

-S->R  works in all angular mode
-But in DEG mode, the longitude L & the latitude B must be expressed in decimal degrees
 

Rectangular->Spherical Conversion
 

-This program employs the same formulae as above

       X = R cos B cos L
       Y = R cos B sin L
       Z = R sin B
 
 

STACK	INPUTS	OUTPUTS
Level 3	X	R
Level 2	Y	L
Level 1	Z	B

Example:             X = 3   Y = 5   Z = 7

       3  ENTER
       5  ENTER
       7  R->S      >>>>   R = 9.11043357915
                                      L = 59°0362434679                ( in DEG mode )
                                      B = 50°2059345082

Notes:

-R->S  works in all angular mode
-But in DEG mode, the longitude L & the latitude B are expressed in decimal degrees
 

Equatorial<->Ecliptic Conversions
 

Formula:

                       sin  b = cos e  sin d - sin e  cos d  sin a
             cos b  cos l = cos d  cos a
             cos b  sin  l = sin e  sin d  +  cos e  cos d  sin a
 
 
 

STACK	INPUTS	OUTPUTS
Level 3	A	/
Level 2	D	L
Level 1	e	B

    Where  A = right-ascension , D = declination , e = obliquity of the ecliptic

Example:      A = 123°    D = 41°   e = 23°439279

          123        ENTER
           41         ENTER
   23.439279       EE          >>>>    L = 116°028898378
                                                      B = 20°4963945697

Notes:

-EE  works in all angular mode
-But in DEG mode, all the angles must be expressed in decimal degrees

-To get the reverse conversion, place -e in level 1
 

AD -> LB
 

 AD->LB   takes the data in TADXYZ and stores in 'TLBXYZ'  the same data
 after performing the equatorial-ecliptic conversion
 

Note:

-This routine is automatically called by GAUSSHG , HERGET & IMPROVE.
 

Warning:

-If flag 1 is clear, the obliquity of the ecliptic at the epoch J2000.0 is used.
-If flag 1 is set, this program uses the obliquity of the ecliptic of the date.
 

Obliquity of the Ecliptic
 

 OBL takes the number of days since 2000/01/01 12h  TT and returns the mean obliquity of the ecliptic in decimal degrees.
 
 

STACK	INPUT	OUTPUT
Level 1	days	obl

Example:      days = 5512.25 = 2015/02/03  18h

     5512.25   OBL  >>>>   obl = 23°4373159804
 

Polar->Rectangular Conversion
 

Formulae:

       X = R cos L
       Y = R sin L
 
 

STACK	INPUTS	OUTPUTS
Level 2	R	X
Level 1	L	Y

Example:      R = 12   L = 128°71

      12       ENTER
  128.71     R->P          X= -7.50454629137                ( in DEG mode )
                                    Y =  9.36385524027

Notes:

-R->P  works in all angular mode
-But in DEG mode, the angles must be expressed in decimal degrees
 

Rectangular->Polar Conversion
 

Formulae:

       X = R cos L
       Y = R sin L
 
 

STACK	INPUTS	OUTPUTS
Level 2	X	R
Level 1	Y	L

Example:        X= -7.50454629137      Y =  9.36385524027

  -7.50454629137   ENTER
   9.36385524027    P->R      >>>>         R = 12
                                                                L = 128°71                ( in DEG mode )

Notes:

-P->R  works in all angular mode
-But in DEG mode, the angles must be expressed in decimal degrees
 
 
 

References:

[1]  "SOLEX" which may be downloaded from http://chemistry.unina.it/~alvitagl/solex/
[2]  "FUNDAMENTALS OF ASTRODYNAMICS" - Roger R. Bate Donald D. Mueller Jerry E. White
[3]  "Fundamentals of Celestial Mechanics" - JMA Danby - Willmann-Bell - ISBN 0-943396-20-4
[3]  Robin M. Green - "Spherical Astronomy" - Cambridge University Press - ISBN  0-521-31779-7
[4]  John D. Anderson - JPL Technical Report n° 32-497
[5]  Dan Boulet - "Methods of Orbit Determination for the Micro Computer" - Willmann-Bell  -  ISBN 0-943396-34-4