Orbit Determination for the HP-41

Overview
 

 0°)  3 Subroutines: "R-S"  "S-R"  "EE"
 1°)  Gauss-Herrick-Gibbs Method
 2°)  Position & Velocity >>> Orbital Elements
 3°)  Position & Velocity >>> Position2
 4°)  2 Positions >>> Velocity
 5°)  Herget's Method: more than 3 Observations
 6°)  A Test
 
 

0°)  3 Subroutines
 

-These short routines are already listed in "Transformation of Coordinates for the HP-41"
 

    a) 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  = ( x² + y² + z² )^1/2  ,   l  = longitude ,  b = latitude

-However, the results can be obtained more easily by the R-P and P-R functions:
 T-register is saved and no data register is used!
 
 

 
( 16 bytes / SIZE 000 )
 
 

 
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     l  = 53.13010235°
                  RDN     b = -54.46232221°
 

   b) Spherical-Rectangular conversion:
 
 

 
( 16 bytes / SIZE 000 )
 
 

 
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

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

   c) A very useful subroutine: "EE"
 

-Many transformations use the same type of formulae which appear in the equatorial-ecliptic conversion, namely:

                       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

- But once again, P-R and R-P lead to a shorter and faster algorithm.
 
 

 
( 31 bytes / SIZE 000 )
 

-If "EE" is useful for you but "R-S" and "S-R" are not, here is another version of this program that doesn't need any subroutine:
 
 

 
( 32 bytes / SIZE 000 )
 
 

  where  e = obliquity of the ecliptic

Example:     if right-ascension = RA = 116.328942 , declination = decl = 28.026183 and  e = 23.4392911

    23.4392911  ENTER^
    28.026183    ENTER^
  116.328942    XEQ "EE"  >>>> l  =  113.215630°   RDN   b  =  6.684170°
 

-Like "R-S" and "S-R" , "EE"  works in all angular modes.
 

1°)  Gauss-Herrick-Gibbs Method
 

-This program is very similar to the routine listed in "Comets for the HP-41" paragraph  6-b)
-However, it does not use any M-Code function, it also works if e = 1 ( parabolic orbits )
 and the velocity vector at the second instant V₂ is computed too.

-Moreover, Newton's method is used instead of the secant method for faster convergence.

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

Data Registers:              R00 & R19 to R50: temp               ( Registers R01 thru R18 are to be initialized before executing "GAUSS" )

                                      •  R01 = t₁               •  R07 = t₂               •  R13 = t₃                  expressed in days since 2000/01/01  0h TT
                                      •  R02 = RA₁           •  R08 = RA₂          •  R14 = RA₃              in hh.mnss
                                      •  R03 = d₁              •  R09 = d₂              •  R15 = d₃                 in °. ' ''
                                      •  R04 = X₁             •  R10 = X₂             •  R16 = X₃                 in Astronomical Units          rectangular ecliptic
                                      •  R05 = Y₁             •  R11 = Y₂             •  R17 = Y₃                 ---------------------              coordinates
                                      •  R06 = Z₁             •  R12 = Z₂              •  R18 = Z₃                 ---------------------               of the Sun

    >>>> When the program stops,

                 R19 = T = instant of passage in perihelion  ( in days from 2000/01/01  0h )
                 R20 = q = perihelion distance    ( in Astronomical Units )
                 R21 = e = eccentricity
                 R22 =  i = inclination  ( in degrees )
                 R23 = omega = argument of perihelion ( in degrees )
                 R24 = OMEGA = longitude of the ascending node ( in degrees )
                 R25 = n = mean motion  ( in degrees/day )
                 R26 = p = parameter ( in AU )

   >>>>  You have also:       R00 = Distance Earth-Comet at the 2nd instant

                 R29 = Distance Sun-Comet at the 1st instant  = r₁
                 R28 = (Distance Sun-Comet at the 2nd instant)^3 = r₂^3
                 R30 = Distance Sun-Comet at the 3rd instant      = r₃

      The coordinates of the vector r₁ are stored in R41-R42-R43
      The coordinates of the vector r₂ are stored in R44-R45-R46         The coordinates of the vector V₂ are stored in R34-R35-R36
      The coordinates of the vector r₃ are stored in R47-R48-R49

Flags:  F00 - F01       Set   F01 if you are using the mean ecliptic of the date(s)
                                  Clear  F01 if you are using the standard equinox J2000
Subroutines: /
 
 
 
 

 
          ( 1102 bytes / SIZE 051 )
 
 

  Where  D  is an estimtion of the distance Sun-Comet at the 2nd 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                 5510                                   5519                                5529                               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
     Y        -0.736974521                     -0.624900515                 -0.482202751                   Astronomical Units
     Z          0.000019390                      0.000019018                   0.000014378                   Astronomical Units
 

-Store these 18 numbers into

                        R01                                     R07                                   R13
                        R02                                     R08                                   R14
                        R03                                     R09                                   R15                               respectively
                        R04                                     R10                                   R16
                        R05                                     R11                                   R17
                        R06                                     R12                                   R18

-If you choose  D = 3 AU

  CF 01   3   XEQ "GAUSS"   >>>>    T = 5536.68812 = R19                                 ---Execution time = 58s---

 and               q      =   2.002314                    = R20
                      e      =   0.999456                    = R21
                       i      = 105°2133                      = R22
                 omega   = -67°7290                      = R23
               OMEGA = -29°5133                      = R24
                      n      =   0.000004415 deg/d     = R25
                      p      =  4.003539                     = R26
 

Notes:

-If your guess leads to a negative value for the distance Sun-asteroid or Earth-asteroid,
 the HP-41 emits a TONE 9 , displays  D=?  and stops line 378
-Key in another guess in X-register and R/S

-If you add  STO 33 after line 422 and if you add  STO 32  after line 402,
 you will get the distance Earth-Asteroid at the 1st & 3rd instants in R32 & R33

-This will be useful in the Herget method below and for aberration too...

-Lines 523 to 727 may be replaced by    RCL 07   XEQ "RV-EL"
  where  "RV-EL" is listed in paragraph 3 below.
-The results could be slightly different because "GAUSS" employs
  the positions at the first and third instant.

-To get accurate position of the Sun, see reference [1] which uses DE421 ephemerides
-If the coordinates of the asteroid are topocentric, use the topocentric coordinates of the Sun.
 

2°)  Position & Velocity >>> Orbital Elements
 

-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
-The formulas are given in "Comets for the HP-41"
 

Data Registers:              R00 = t                    ( Registers R44-45-46 & R34-35-36 are to be initialized before executing "RV-EL" )

                                      •  R44-R45-R46 = r  = rectangular ecliptic coordinates of the asteroid
                                      •  R34-R35-R36 = V = -------------------------------------  velocity    at the same instant t

    >>>> When the program stops,

                 R19 = T   ( in days from 2000/01/01  0h )
                 R20 = q    ( in Astronomical Units )
                 R21 = e
                 R22 =  i     ( in degrees )
                 R23 = omega  ( in degrees )
                 R24 = OMEGA  ( in degrees )
                 R25 = n = mean motion  ( in degrees/day )
                 R26 = p = parameter
                 R27 = a = semi-major axis
 

Flag:  F00
Subroutines: /

-Line 100  is a three-byte GTO 02
 
 

 
        ( 363 bytes / SIZE 049 )
 
 

  Where  t  &  T  are expressed in days since  2000/01/01  0h  TT

Example1:                t = 4321  r = ( 1.4  ,  5.3  ,  -0.9 ) V = ( 0.003 , -0.004 , -0.009 )

   Store  the 3 components of r into  R44-R45-R46
     and   -------------------- V --   R34-R35-R36   respectively
 

   4321   XEQ "RV-EL"  >>>>   T = 4487.011977 = R19                   ---Execution time = 12s---

  and in registers R20 thru R27

                      q      =   5.419995                    = R20
                      e      =   0.990189                    = R21
                       i      = 112°36768                    = R22
                 omega   = -151°91629                  = R23
               OMEGA = -100°92280                  = R24
                      n      =   0.000075905 deg/d     = R25
                      p     =   10.786814                   = R26
                      a      =  552.446418                 = R27

Example2:       Same values except the last component of the velocity:      -0.010    STO 36

   4321   XEQ "RV-EL"  >>>>   T = 4431.72425 = R19                    ---Execution time = 12s---

  and in registers R20 thru R27

                      q      =   5.474724                    = R20
                      e      =   1.341612                    = R21
                       i      = 110°43073                    = R22
                 omega   = -157°13432                  = R23
               OMEGA = -101°29046                  = R24
                      n      =   -0.015362 deg/d         = R25
                      p     =   12.819681                   = R26
                      a      =  -16.026128                 = R27

Notes:

-For a hyperbolic orbit,  n & a are negative
-For a parabolic orbit, n = 0
 

3°)  Position & Velocity >>> Position2
 

-Here, we could use "RV-EL" first to get the orbital elements and then the programs listed in "Comets for the HP-41" to get the position at another instant.
-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
 

Data Registers:              R00 = t₂ - t₁                    ( Registers R44-45-46 & R34-35-36 are to be initialized before executing "RV-R" )

                                      •  R44-R45-R46 = r₁  = rectangular ecliptic coordinates of the asteroid
                                      •  R34-R35-R36 = V₁ = -------------------------------------  velocity    at the same instant t₁

                                          R01 thru R06 & R12 thru R26 are unused

    >>>> When the program stops,

                 R27 = X    rectangular  ecliptic
                 R28 = Y    heliocentric coordinates                                                 R27-R28-R29 = r₂
                 R29 = Z     of the asteroid at the second instant in AU.

Flag:  F00   ( F00 is temporarily set if the orbit is hyperbolic )
Subroutines: /
 

-Line 183 is a three-byte  GTO 01
 
 

 
     ( 310 bytes / SIZE 047 )
 
 

  Where  t₂ - t₁   is expressed in days and  X Y Z  are the rectangular ecliptic coordinates of the asteroid ( in AU )

Example1:        r₁  = ( 0.16 , 1.38 , 0.24 )    V₁ = ( 0.015 , 0.01 , 0.001 )    Find the position 100 days later

     0.16   STO 44     0.015   STO 34
     1.38   STO 45     0.010   STO 35
     0.24   STO 46     0.001   STO 36

     100   XEQ "RV-R"   >>>>    X = 1.509299637 = R27                              ---Execution time = 25s---
                                      RDN     Y = 1.919542031 = R28
                                      RDN     Z = 0.265117223 = R29

Example2:      Same values except the second component of the velocity vector:  0.015 instead of 0.01

    0.015   STO 35
     100      R/S       >>>>   X = 1.541288717 = R27                              ---Execution time = 31s---
                             RDN     Y = 2.468789822 = R28
                             RDN     Z = 0.277102516 = R29
 

Notes:

-The HP-41 displays the successive approximations
-In example 1 , the orbit is elliptic.
-In example 2 , the orbit is hyperbolic.

-I first wrote a program that uses the initial guess given in reference [2]
-It works well for elliptic orbits, but for hyperbolic orbits, there are a few iterations only if x is large.
-When "RV-R" is called by "HRG" below, I found better to take as a 1st guess x₀ = k ( t₂ - t₁ ) / (-a)
-If you prefer the other value, replace lines 44 to 51 by
 
 

 
-Even simpler:  delete lines 45 to 51 and add  SIGN  STO 10 after line 03
 

4°)  2 Positions >>> Velocity
 

-We solve Lambert problem: given 2 position-vectors  r₁ & r₂  at different instants  t₁ & t₂
  "RR-V"  returns the velocity vector  V₁ at the first instant
-It uses an iteration method called the method of Herrick & Liu  or  p-iteration method ( cf reference [3] )

-You have to give an estimate of the parameter p

-Flag F04 must be cleared if the angle  ( r₁ ; r₂ ) < 180°
-Flag F04 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

-This equation is solved by the secant method.
 
 

Data Registers:                                             ( Registers R44-45-46 & R41-42-43 are to be initialized before executing "RV-R" )

                                      •  R44-R45-R46 = r₁  = rectangular ecliptic coordinates of the asteroid at the instant t₁
                                      •  R41-R42-R43 = r₂ = ----------------------------------------------- at the instant t₂

    >>>> When the program stops,   R34-R35-R36 = V₁ = rectangular ecliptic coordinates of the velocity at the first instant t₁
 

Flag:  F04    Clear F04 if the angle  ( r₁ ; r₂ ) < 180°
                     Set F04 if the angle  ( r₁ ; r₂ ) > 180°

Subroutine:  "RV-R"  ( cf paragraph 3 above )
 
 
 

 
      ( 286 bytes / SIZE 049 )
 
 

   When the program stops,  R03 = p

Example:    t₁ = 0   ,   t₁ = 100 days  ,    r₁ = ( 0.16 , 1.38 , 0.24 )  ,  r₂ = ( 1.509299637 , 1.919542031 , 0.265117223 )

     0.16   STO 44     1.509299637   STO 41
     1.38   STO 45     1.919542031   STO 42
     0.24   STO 46     0.265117223   STO 43

-If you choose  p₀ = 1   and   CF 04

           CF 04
     0    ENTER^
   100  ENTER^
     1    XEQ "RR-V"  >>>>   p = 1.276338017                                         ---Execution time = 2mn58s---

 and the velocity-vector in R34-R35-R36 = ( 0.015 , 0.01 , 0.001 )
 

Notes:

-The HP-41 displays the successive p-approximations and, between 2 such approximations, the successive x-approximations !

-A bad guess may lead to a DATA ERROR line 112.
-Try another p-estimation ( often a smaller one ) and execute "RR-V" again.

-After finding the velocity-vector, you can use "RV-EL" to get the orbital elements...
-You will find for instance eccentricity = e = 0.771671...
-The parameter p in R26 will be sometimes slightly different than the value found above:
-In this example R26 = 1.276338002 instead of  1.276338017
 

5°)  Herget's Method: more than 3 Observations
 

-Now, we have the tools to program Herget's method with more than 3 observations:

-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-41 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
 
 

Data Registers:              R00 & R19 to R50: temp               ( Registers R50 thru R50+6n are to be initialized before executing "HRG" )

                                      •  R50 = control number of the data below = 51.00006 + ( 50+6n) / 1000

                                      •  R51 = t₁               •  R57 = t₂         ...............    •  R45+6n = t_n                  expressed in days since 2000/01/01  0h TT
                                      •  R52 = L₁              •  R58 = L₂        ...............    •  R46+6n = L_n                Ecliptic longitudes in degrees
                                      •  R53 = b₁              •  R59 = b₂         ...............   •  R47+6n = b_n                 Ecliptic latitudes in degrees
                                      •  R54 = X₁             •  R60 = X₂         ..............    •  R48+6n = X_n                 in Astronomical Units          rectangular ecliptic
                                      •  R55 = Y₁             •  R61 = Y₂         ..............    •  R49+6n = Y_n                 ---------------------              coordinates
                                      •  R56 = Z₁             •  R62 = Z₂          ..............    •  R50+6n = Z_n                 ---------------------               of the Sun

    >>>> When the program eventually stops,

                 R19 = T = instant of passage in perihelion  ( in days from 2000/01/01  0h )
                 R20 = q = perihelion distance    ( in Astronomical Units )
                 R21 = e = eccentricity
                 R22 =  i = inclination  ( in degrees )
                 R23 = omega = argument of perihelion ( in degrees )
                 R24 = OMEGA = longitude of the ascending node ( in degrees )
                 R25 = n = mean motion  ( in degrees/day )
                 R26 = p = parameter ( in AU )
 

Flags:  F00 - F04      Clear F04 if the angle  ( r₁ ; r_n ) < 180°
                                     Set F04 if the angle  ( r₁ ; r_n ) > 180°

Subroutines:  "RR-V"  "RV-R"  "S-R"
 
 

-Line 148 is a three-byte  GTO 10
 
 

 
      ( 471 bytes / SIZE 69+6n )
 
 

 
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                 5510                         5519                         5529                     5538                     5547                       5557
    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 these 36 numbers into

                        R51                        R57                            R63                                                                                    R81
                        R52                        R58                            R64                                                                                    R82
                        R53                        R59                            R65                ..................................................                  R83
                        R54                        R60                            R66                                                                                    R84
                        R55                        R61                            R67                                                                                    R85
                        R56                        R62                            R68                                                                                    R86

-Before using "HRG" , let's execute "GAUSS" with the first 3 observations and with the last 3 observations

  51.001018  REGMOVE  CF 01
   3  ( or another guess )  XEQ "GAUSS"   >>>>      T = 5536.50275 = R19                                 ---Execution time = 58s---

 and               q      =   2.004237                    = R20
                      e      =   1.004567                    = R21
                       i      = 105°2001                      = R22
                 omega   = -67°7889                      = R23
               OMEGA = -29°4726                      = R24
                      n      =  -0.000107203 deg/d     = R25
                      p      =  4.017627                     = R26     and   R00 = 2.405567

and, if you have modified "GAUSS" as suggested in the 1st paragraph

 R32 = D₁ = 2.316488230

-Likewise, for the last 3 observations

  69.001018  REGMOVE
   3  ( or another guess )  XEQ "GAUSS"   >>>>      T = 5536.66386 = R19                                 ---Execution time = 58s---

 and               q      =   2.004552                    = R20
                      e      =   1.003666                    = R21
                       i      = 105°1992                      = R22
                 omega   = -67°7051                      = R23
               OMEGA = -29°4677                      = R24
                      n      =  -0.000077070 deg/d     = R25
                      p      =  4.016452                      = R26     and   R00 = 2.655687

and, if you have modified "GAUSS" as suggested in the 1st paragraph

 R33 = D_n = 2.717138781

-Note that we get a hyperbolic orbit, unlike the result found in paragraph 1 without rounding the data.

2nd Step         Now, we have to replace the right-ascensions & declinations by the longitudes & latitudes.
                       This may be performed by the short routine below:
 
 

 
     ( 101 bytes )
 
 

     where n = number of observations

Example:     Here  n = 6  so,

  CF 01    6  XEQ "INIT"   >>>>   51.08606 = R50

[ CF 01 if the coordinates are reffered to the standard equinox  J2000
  SF 01 ------------------------------------  equinox of the date ]
 

3rd Step     Now, we are ready to use "HRG"

-A good choice for h is often  0.01  AU
-The results given by "GAUSS"   suggest  p = 4.017   D₁ = 2.316488230   D_n = 2.717138781  thus,

                           CF 04

         0.01           ENTER^
        4.017          ENTER^
  2.316488230    ENTER^
  2.717138781    XEQ "HRG"  >>>>    D₁ = 2.314988382 = R05              ---Execution time = 10m55s---
                                                X<>Y    D_n = 2.715024458 = R06

-To get the next approximation, press R/S  without  clearing X,  it yields

                           R/S    >>>>   D₁ = 2.314978136 = R05
                                    X<>Y   D_n = 2.715012525 = R06

                           R/S    >>>>   D₁ = 2.314977837 = R05
                                    X<>Y   D_n = 2.715012172 = R06

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

4th Step       If we want to stop the iterations and get the orbital elements, clear register X and press R/S

             CLX   R/S   >>>>          T = 5536.68133 = R19

 and               q      =   2.002584                    = R20
                      e      =   1.000091                    = R21
                       i      = 105°21130                    = R22
                 omega   = -67°7278                      = R23 = 292°2722
               OMEGA = -29°5072                      = R24 = 330°4928

-The exact values - given for instance by the imcce - are

              T =  2015 février 27,64787 TT ±0,00000 TT = 5536.64787  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
 

Notes:

-Due to the execution time, a good emulator or a 41CL is preferable...
-The results seem better than what we got with "GAUSS" in paragraph1 with more accurate data.
-Convergence would have occured even with the estimates  p = D₁ = D_n =  2  but good guesses will reduce the number of iterations.
 

6°)  A Test
 

-This routine calculates the sum of the squares of the angular separations between the observed positions and the calculated positions.

-We assume that registers R44-R45-R46 = Position-vector & R34-R35-R36 = Velocity-vector
-The n observations data are to be stored in registers R51-R52- ..............

-The right-ascensions & declinations must be replaced by the ecliptic longitudes & latitudes
  ( Execute "INIT" listed above if it's not already done )
 

Data Registers:              R00 & R19 to R50: temp               ( Registers R50 thru R50+6n are to be initialized before executing "EPS" )

                                      •  R44-R45-R46 = r  = rectangular ecliptic coordinates of the asteroid
                                      •  R34-R35-R36 = V = -------------------------------------  velocity    at the same instant t

                                      •  R50 = control number of the data below = 51.00006 + ( 50+6n) / 1000

                                      •  R51 = t₁               •  R57 = t₂         ...............    •  R45+6n = t_n                  expressed in days since 2000/01/01  0h TT
                                      •  R52 = L₁              •  R58 = L₂        ...............    •  R46+6n = L_n                Ecliptic longitudes in degrees
                                      •  R53 = b₁              •  R59 = b₂         ...............   •  R47+6n = b_n                 Ecliptic latitudes in degrees
                                      •  R54 = X₁             •  R60 = X₂         ..............    •  R48+6n = X_n                 in Astronomical Units          rectangular ecliptic
                                      •  R55 = Y₁             •  R61 = Y₂         ..............    •  R49+6n = Y_n                 ---------------------              coordinates
                                      •  R56 = Z₁             •  R62 = Z₂          ..............    •  R50+6n = Z_n                 ---------------------               of the Sun
 

Flags:  F02      Clear F02 if the instant t is the 1st one, i-e  in R51
                          Set   F02 ------------------  2nd one, i-e in R57

Subroutines:   "RV-R"  "R-S"
 
 

 
     ( 100 bytes / SIZE 51+6n )
 
 

     where  eps = SUM [ dL²Cos²b + db² ]

Example:   After solving the example of the previous paragraph ( "HRG" ),

      CF 02  XEQ "EPS"   >>>>   eps = 34  E-9

-Since there are 6 observations, the Root-Mean-Square error is obtained by

     6  /   SQRT   which yields    RMS = 0.000075 deg   about  0.27 arcsecond

Notes:

-If the position vector & the velocity vector were obtained by "GAUSS" and the first 3 observations, set flag F02
 ( don't forget to execute "INIT" before "EPS" )
-If you want to get directly the RMS error, add the following instructions after line 49:

    RCL 50   FRC  ISG X   INT  /   SQRT

-The formula  dL²Cos²b + db² for the square of the angular separations is only valid if the separation is small.
 

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