Comets for the HP-41

Overview
 

 1°)  Olbers' Method  ( Determination of a Parabolic Orbit )
 2°)  Parabolic Motion
 3°)  Elliptic Motion
 4°)  Hyperbolic Motion
 5°)  Geocentric Coordinates ( Subroutine )
 6°)  Gauss' Method  ( Determination of Elliptic, Parabolic and Hyperbolic Orbits )

   a)  Program#1
   b)  Program#2
 

-In the following programs, we have used these notations:

       T = time of passage in perihelion                 i      =  inclination
       q = perihelion distance                          omega   = argument of the perihelion
       e = eccentricity                                   OMEGA = longitude of the ascending node
 

1°) Olbers' Method  ( Determination of a Parabolic Orbit )
 
 

-Olbers' method computes approximate values of  T , q , i , omega , OMEGA  from 3 observations of a comet.
-We assume that the orbit is - at least nearly - parabolic:  e = 1
-The method fails if the inclination i = 0 or is very small.

-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₁
-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₂
-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₃

-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 should be 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

-The ecliptic longitudes and latitudes of the comet are  L_i = Atan m_i/l_i  ( use R-P to get the correct quadrant ) and  b_i = Asin n_i
-Let D_i = distance Earth-Comet at the instant t_i , the D_i 's are unknown but   Q = D₃/D₁ may be calculated by the approximate formula:

    Q = ((t₃ - t₂)/(t₂ - t₁)) ( sin b₁ / sin b₃ ) { [ sin ( L₂ - L'₂ ) / tan b₂ - sin ( L₁ - L'₂ ) / tan b₁ ] / [ sin ( L₃ - L'₂ ) / tan b₃ - sin ( L₂ - L'₂ ) / tan b₂ ] }

    where  L'₂  is the ecliptic longitude of the Sun at the instant  t₂

-"OLBERS" uses the equivalent formula:   Q = - ((t₃ - t₂)/(t₂ - t₁)) [ u₁. ( u₂ x R₂ ) ] / [ u₃. ( u₂ x R₂ ) ]       where  u_i ( l_i , m_i , n_i )  and  R₂ ( X₂ , Y₂ , 0 )

-The rectangular heliocentric coordinates of the comet at the instant   t₁ &  t₃  are then:

      x₁ = l₁ D₁ - X₁               x₃ = Q l₃ D₁ - X₃
      y₁ = m₁ D₁ - Y₁             y₃ = Q m₃ D₁ - Y₃
      z₁ = n₁ D₁                      z₃ = Q n₃ D₁

-So the square of the distance  c² = (x₃ - x₁)² + (y₃ - y₁)² + (z₃ - z₁)²    may be written under the form:   c² = A D₁² + B D₁ + C

-On the other hand:

    r₁² = D₁² - 2 D₁ ( l₁ X₁ + m₁ Y₁ ) + X₁² + Y₁²                             where  r_i = radius vector of the comet at the instant  t_i
    r₃² = Q² D₁² - 2.Q D₁ ( l₃ X₃ + m₃ Y₃ ) + X₃² + Y₃²

-We have also   c = ( r₁ + r₃ ) sin 2µ

   with  µ = Asin { sqrt(2) Sin [ (1/3) Asin [ 3k ( t₃ - t₁ ) ( r₁ + r₃ ) ^-3/2 /sqrt(2) ] ] }     where  k = 0.01720209895 = Gaussian gravitational constant

>>> So  D₁ is the solution of a relatively complicated equation that the following program solves by the secant method.

-When D₁ is found, we compute the quantities

    s₁ = x₁ y₃ - x₃ y₁
    s₂ = y₁ z₃ - y₃ z₁          and    S = ( s₁² + s₂² + s₃² )^1/2
    s₃ = x₁ z₃ - x₃ z₁

-The longitude of the ascending node is then   OMEGA = Atan2 s₂/s₃  ( use R-P to get the proper quadrant )
                   and the inclination                               i       = Acos s₁/S

-    ( v_i + omega ) may be obtained by

    r_i Sin ( v_i + omega ) = z_i / Sin i
    r_i Cos ( v_i + omega ) = x_i Cos OMEGA + y_i Sin OMEGA

-The relation  v₃ - v₁ = Asin S/(r₁ r₃)  is also used - we assume that  v₃ - v₁   does not exceed  90°

-The perihelion distance  q = [ r₁ r₃ Sin² (v₃ - v₁)/2 ] / [ r₁ + r₃ - 2(r₁ r₃)^1/2 Cos (v₃ - v₁)/2 ]

-  (v₃ + v₁) is calculated with the R-P function from the relations

   2.S Sin (v₃ + v₁)/2 = q ( r₃ - r₁ ) Cos (v₃ - v₁)/2
   2.S Cos (v₃ + v₁)/2 = [ q ( r₃ + r₁ ) - r₁ r₃ ] Sin (v₃ - v₁)/2

   whence  v₁  whence  omega since we already know ( v_i + omega )

-Finally    T = t₁ - sqrt(2)/(3k) q^3/2 ( 3.s + s³ )    with  s = Tan v₁/2
 
 

Data Registers:              R00 & R16 to R29: temp        ( Registers R01 thru R15 are to be initialized before executing "OLBERS" )

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

    >>>> When the program stops,

                                        R16 = T   ( in days from 2000/01/01  0h )
                                        R17 = q    ( in Astronomical Units )
                                        R18 =  i     ( in degrees )
                                        R19 = omega  ( in degrees )
                                        R20 = OMEGA  ( in degrees )
Flag:  F00
Subroutines: /
 

-The starting value is 1 AU ( line 66 )
-If the process seems to diverge, stop the program, set flag F00
  and store another guess - for instance 3 - into register R00
-Then, press  XEQ 16
-Line 152  this constant = sqrt(2)/(3k)   where  k = 0.01720209895 = Gaussian gravitational constant
 
 

 
    ( 562 bytes / SIZE 030 )
 
 

 
Example:     You observe a comet during November 2007 and you measure its position:

        2007/11/21  0h TT                  2007/11/24   0h  TT                  2007/11/27   0h  TT
          ( t = 2881 days )                       ( t = 2884 days )                       ( t = 2887 days )                   since  2000/01/01 0h  TT

        RA = 17h07m25s03                RA = 17h07m11s94                  RA = 17h06m59s04
       Decl = -34°21'46"09               Decl = -35°53'19"78                 Decl = -37°24'57"39

-If you use "SUN2" ( cf "Astronomical Ephemeris for the HP-41" §8 )  and P-R , you get the rectangular ecliptic coordinates of the Sun:

         X = -0.519356                        X = -0.473907                           X = -0.427153
         Y = -0.840463                        Y = -0.866219                           Y = -0.889590

-Store these 15 numbers into

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

  and  XEQ "OLBERS"   the successive  D₁-values are displayed ( the last displayed value is 1.867064 ) and eventually:

                      T      =  2902.471673                = R16                   --- execution time = 91 seconds ---
    RDN          q      =   0.969357                     = R17
    RDN           i      =   117°6681                     = R18
    RDN     omega   = -126°3586 = 233°6414 = R19
 LASTX  OMEGA =  111°5459                     = R20

-If we reduce these elements to the standard epoch J2000 ( cf for instance "Orbital Elements from one Equinox to another one" ) it yields:

                      i₂₀₀₀        = 117°6686
                 omega₂₀₀₀   = 233°6403
               OMEGA₂₀₀₀ = 111°4352

-This comet is actually  C/2007 T1 mcNaught  and the exact elements are:
 

            T       = 2902.49731 = 2007/12/12.49731
            q        =  0.969480
            e        =  1.000785                    the orbit is a hyperbola
            i         = 117°6490
        omega    =  233°6712                   referred to J2000.0
      OMEGA  =  111°4186
 

-The error is about 37 minutes for the time of passage in perihelion, q is good and the angle errors are smaller than 0.03°
-Our results are satisfactory all the more that we have neglected the aberration, parallax and planetary perturbations!
 

2°)  Parabolic Motion
 

-Let  s = Tan v/2  where  v = true anomaly,
 s  is the unique solution of the cubic equation:   s³ + 3s - W = 0
 where  W = (3k/sqrt(2)) (t-T) q ^-3/2   and  k = 0.01720209895 = Gaussian gravitational constant

-Then, the radius vector is obtained by  r = q(1+s²)
 
 

  >>> Registers R01 thru R06 are to be initialized before executing "PARM"
 

Data Registers:   R00 = t ( in days at the beginning , in millenia at the end ) from 2000/01/01  0h TT

                            •  R01 = T ( in days from 2000/01/01  0h TT )       •  R04 = i
                            •  R02 = q                                                             •  R05 = omega
                        (   •  R03 = e = 1    this register is actually unused )     •  R06 = OMEGA

                              R07 = RA ( hh.mnss )                           R10 = l_g = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
                              R08 = decl ( °. ' '' )                               R11 = b_g = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
                              R09 = r_T = distance Earth-Comet ( AU )                                                                                R14 = b = heliocentric latitude ( in degrees )

                              R15 = X  &  R16 = Y  ( rectangular coordinates of the Sun )

                >>>>    All the coordinates are referred to the mean ecliptic and equator of the date.

Flags: /
Subroutines:  "GEO"                         ( cf paragraph 5 )
                         "SUN2" or "SUN3"  ( cf "Astronomical Ephemeris for the HP-41" )
                          "S-R" "R-S" "EE" ( cf "Transformation of Coordinates for the HP-41" )

-Line 06 , this constant is  sqrt(8)/(3k)  where  k =  0.01720209895 = Gaussian gravitational constant
 
 

 
   ( 62 bytes / SIZE 017 )
 
 

  with t expressed in days from 2000/01/01  0h  TT

Example:     Calculate the geocentric position of comet Kohler for 1977 September 29.0 TT using the following elements

    T = 1977 November 10.5659  TT         i₂₀₀₀      =   48°7131
    q =  0.990662                                 omega₂₀₀₀   = 163°4788
    e =  1                                            OMEGA₂₀₀₀ = 182°1660

-First, we reduce the elements to 1977/11/29  it gives:

          i       =    48°7160
     omega   =  163°4793     ( with the programs "IOOM" or "IOO" )
   OMEGA = 181°8571

-We have  T = -8086.4341 days  and  t = -8129 days  from  2000/01/01 0h  TT  so we store

    -8086.4341             48.7160                                         R01            R04
     0.990662              163.4793                into                   R02            R05               respectively
     1                           181.8571                                         R03            R06

-Then:

         -8129  XEQ "PARM"  >>>>   RA = 16h19m09s          -- execution time = 24 seconds ---
                                             RDN   decl =  20°16'25"
                                             RDN     r_T  =  1.3062 AU
                                             RDN    Psi  =  62°51

-If you want to get the coordinates referred to J2000 , use for instance "PREQ" ( cf "Transformation of coordinates for the HP-41" ) , it yields:

        RA₂₀₀₀ =  16h20m08s
       Decl₂₀₀₀ =  20°13'15"

-If you have to calculate several positions, i , omega , OMEGA may be regarded as constant for several days or weeks
 ( though "IOO" is relatively fast to use... )
 

3°)  Elliptic Motion
 

-We have to solve Kepler's equation   M = E - e sin E  ( in radians )     where  M = k(t-T) ((1-e)/q)^3/2  = mean anomaly and  E = excentric anomaly
-This equation is solved by Newton's method.
-Then, the true anomaly v is obtained by   Tan v/2 = ((e+1)/(1-e))^1/2  Tan E/2
  and the radius vector is   r = q(1+e)/(1+e cos v)
 
 

  >>> Registers R01 thru R06 are to be initialized before executing "ELLM"
 

Data Registers:   R00 = t ( in days at the beginning , in millenia at the end ) from 2000/01/01  0h TT

                            •  R01 = T ( in days from 2000/01/01  0h TT )       •  R04 = i
                            •  R02 = q                                                             •  R05 = omega
                            •  R03 = e < 1                                                       •  R06 = OMEGA

                              R07 = RA ( hh.mnss )                           R10 = l_g = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
                              R08 = decl ( °. ' '' )                               R11 = b_g = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
                              R09 = r_T = distance Earth-Comet ( AU )                                                                                R14 = b = heliocentric latitude ( in degrees )

                              R15 = X  &  R16 = Y  ( rectangular coordinates of the Sun )

                >>>>    All the coordinates are referred to the mean ecliptic and equator of the date.

Flags: /
Subroutines:   "GEO"                         ( cf paragraph 5 )
                        "SUN2" or "SUN3"  ( cf "Astronomical Ephemeris for the HP-41" )
                        "S-R" "R-S" "EE" ( cf "Transformation of Coordinates for the HP-41" )

-Line 06 , this constant = 1/k  where  k =  0.01720209895 = Gaussian gravitational constant
-Lines 18 thru 73 use the method given by J Meeus in reference [2]
  to find a good starting value for the iterations, even in very difficult cases like  M = 7°  e = 0.999
-The iterations stops when the difference between 2 consecutive approximations is smaller than  0.000001 degree
-Change line 91 if you want a better - or a lesser - accuracy.
 
 

 
     ( 166 bytes / SIZE 017 )
 
 

  with t expressed in days from 2000/01/01  0h  TT

Example:   Calculate the geocentric position of comet C/2007 K6 for 2007 December 01  0h  TT using the following elements
 

    T = 2007 July 01.47533  TT                 i₂₀₀₀      =  105°063204
    q =  3.432968                                 omega₂₀₀₀   = 337°140230
    e =  0.984585                               OMEGA₂₀₀₀ = 298°075386

-We reduce the elements to 2007/12/01  it gives

            i       =  105°06377
       omega   =  337°13933
     OMEGA = 298°18572

-We have  T = 2738.47533 days  and  t = 2891 days  from  2000/01/01 0h  TT  so we store

    2738.47533           105.06377                                       R01            R04
    3.432968               337.13933              into                   R02            R05               respectively
    0.984585               298.18572                                       R03            R06

-Then    2891  XEQ "ELLM"  >>>>   RA = 19h07m27s          -- execution time = 39 seconds ---
                                                RDN   decl = -15°25'02"
                                                RDN     r_T  =  4.4261 AU
                                                RDN    Psi  =  38°52

-Referred to J2000 , it yields:

        RA₂₀₀₀ =  16h20m08s
       Decl₂₀₀₀ =  20°13'15"
 

4°)  Hyperbolic Motion
 

-For hyperbolic orbits, Kepler's equation becomes:    k(t-T) ((e-1)/q)^3/2  = e Sinh E - E
-Newton's method is again used.
-Then, the true anomaly v is obtained by   Tan v/2 = ((e+1)/(e-1))^1/2  Tanh E/2
  and the radius vector by   r = q(1+e)/(1+e cos v)

-The following program may be improved if you have hyperbolic functions in M-Code...
 

  >>> Registers R01 thru R06 are to be initialized before executing "HYPM"
 

Data Registers:   R00 = t ( in days at the beginning , in millenia at the end ) from 2000/01/01  0h TT

                            •  R01 = T ( in days from 2000/01/01  0h TT )       •  R04 = i
                            •  R02 = q                                                             •  R05 = omega
                            •  R03 = e > 1                                                       •  R06 = OMEGA

                              R07 = RA ( hh.mnss )                           R10 = l_g = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
                              R08 = decl ( °. ' '' )                               R11 = b_g = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
                              R09 = r_T = distance Earth-Comet ( AU )                                                                                R14 = b = heliocentric latitude ( in degrees )

                              R15 = X  &  R16 = Y  ( rectangular coordinates of the Sun )

               >>>>    All the coordinates are referred to the mean ecliptic and equator of the date.

Flags: /
Subroutines:   "GEO"                         ( cf paragraph 5 )
                        "SUN2" or "SUN3"  ( cf "Astronomical Ephemeris for the HP-41" )
                        "S-R" "R-S" "EE" ( cf "Transformation of Coordinates for the HP-41" )

-Line 06 , this constant = 1/k  where  k =  0.01720209895 = Gaussian gravitational constant
-The iterations stops when the difference between 2 consecutive approximations is smaller than  0.000001 degree
-Change line 55 if you want a better - or a lesser - accuracy.
 
 

 
   ( 126 bytes / SIZE 017 )
 
 

  with t expressed in days from 2000/01/01  0h  TT

Example:  Calculate the geocentric position of comet C/2007 T1 for 2008 January 01  6h  TT using the following elements
 

    T = 2007 December 12.49731  TT                i₂₀₀₀      = 117°649041
    q =  0.969480                                          omega₂₀₀₀   = 233°671201
    e =  1.000785                                        OMEGA₂₀₀₀ = 111°418623

                                                                            i       =  117°64857
-Reduce the elements to 2008/01/01.25          omega   =  233°67226
                                                                     OMEGA = 111°53088

-We have  T = 2902.49731 days  and  t = 2922.25 days  from  2000/01/01 0h  TT  so we store

    2902.49731           117.64857                                       R01            R04
    0.969480               233.67226              into                   R02            R05               respectively
    1.000785               111.53088                                       R03            R06

-Then    2922.25  XEQ "HYPM"  >>>>   RA = 17h02m54s          -- execution time = 30 seconds ---         Referred to J2000 , it yields:
                                                    RDN   decl = -57°40'49"
                                                    RDN     r_T  =  1.5825 AU                                                                         RA₂₀₀₀ =  17h02m13s
                                                    RDN    Psi  =  39°16                                                                                Decl₂₀₀₀ =  -57°40'09"

>>>> Let's use the approximate elements that we've found in paragraph 1 after reducing  i , omega and OMEGA to the epoch J2008

  XEQ "PARM"  gives   RA = 17h02m52s  ,  Decl = -57°40'13" ,  r_T  =  1.5825 AU ,  Psi = 39°15

-The difference is about 36 arcseconds only, this accuracy is sufficient for telescope pointing.
 

5°)  Geocentric Coordinates
 

-"GEO" is a subroutine called by the programs listed in paragraphs 2-3-4
-It takes the radius vector and the true anomaly in registers X and Y respectively
 and returns the right ascension, declination and distance to the Earth in registers X , Y , Z
-The heliocentric ecliptic coordinates are also computed and stored in registers R12  R13  R14
 

Data Registers:   R00 = t  ( in days at the beginning , in millenia at the end )  from 2000/01/01  0h TT

                              R01 = T ( in days from 2000/01/01  0h TT )         R04 = i
                              R02 = q                                                               R05 = omega
                              R03 = e                                                               R06 = OMEGA

                              R07 = RA ( hh.mnss )                           R10 = l_g = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
                              R08 = decl ( °. ' '' )                               R11 = b_g = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
                              R09 = r_T = distance Earth-Comet ( AU )                                                                                R14 = b = heliocentric latitude ( in degrees )

                              R15 = X  &  R16 = Y  ( rectangular ecliptic coordinates of the Sun )

                >>>>    All the coordinates are referred to the mean ecliptic and equator of the date.

Flags: /
Subroutines:  "SUN2" or "SUN3"  ( cf "Astronomical Ephemeris for the HP-41" )
                        "S-R" "R-S" "EE"      ( cf "Transformation of Coordinates for the HP-41" )

-If you don't want to save the coordinates of the Sun in R15 & R16, replace lines 33-35-40-43 by
  STO 08   STO 09   RCL 09   RCL 08   respectively:
-2 bytes will be saved and SIZE 015 will be enough.
 
 

 
    ( 124 bytes / SIZE 017 )
 
 

      r = radius vector       RA = right-ascension     r_T = distance Earth-Comet             ( Execution time = 21 seconds )
      v = true anomaly       Decl = declination        Psi =  elongation = angle Sun-Earth-Comet

-Examples are given in paragraphs 2-3-4.
 

6°) Gauss' Method  ( Determination of Elliptic, Parabolic and Hyperbolic Orbits )
 

     a)  Program#1
 

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

-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₁
-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₂
-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₃

-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 should be 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₂³   ( 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₃ )                               ( all the Z_i's ~ 0 if we use the rectangular ecliptic coordinates of the Sun )

         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₂³   with    r₂ =  D₂² - 2 D₂ ( l₂ X₂ + m₂ Y₂ ) + X₂² + Y₂²

-Then,  D₁ = ( P₁ + Q₁ / r₂³ )/( A₁ + B₁ / r₂³ )  and  D₃ = ( P₃ + Q₃ / r₂³ )/( A₃ + 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

-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                          you can also use C = r₁ x r₃  instead of  h  to get i and OMEGA

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

-    ( v_i + omega ) is obtained by

    r_i Sin ( v_i + omega ) = z_i / Sin i
    r_i Cos ( v_i + omega ) = x_i Cos OMEGA + y_i Sin OMEGA        whence   (v₃ - v₁)/2

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

    v + omega = sign(z) Acos [ ( x Cos OMEGA + y Sin OMEGA ) / r ]   with perhaps less accurate results in some cases.

-The relations:   2e Sin (v₃ + v₁)/2 = ( p/r₁ - p/r₃ ) / Sin (v₃ - v₁)/2                   give  e and (v₃ + v₁)/2  whence v₁ and omega
                        2e Cos (v₃ + v₁)/2 = ( p/r₁ + p/r₃ - 2 ) / Cos (v₃ - v₁)/2

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

Data Registers:              R00 & R16 to R41: temp                ( Registers R01 thru R15 are to be initialized before executing "GAUSS" )

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

    >>>> When the program stops,

                        R16 = T   ( in days from 2000/01/01  0h )
                        R17 = q    ( in Astronomical Units )
                        R18 = e
                        R19 =  i     ( in degrees )
                        R20 = omega  ( in degrees )
                        R21 = OMEGA  ( in degrees )
                        R22 = n = mean motion  ( in degrees/day ) - virtual for a hyperbolic orbit.

Flag:  F00
Subroutines: /

-Line 110 , this constant = 1/k^2   where  k = 0.01720209895 = Gaussian gravitational constant
-Line 255 , this loop solves an equation by the secant method  ( the unknown is the distance Earth-Comet at the second instant ), starting-value = 1 ( line 217 )
-If the process seems to diverge or if it converges to a negative value:
-Stop the program, set flag F00 , store another guess into R00 and press  XEQ 10.

-Lines 353 to 425 calculate and store the components of the velocity at the second instant into R34-35-36.
-However, these values are only used to compute the angular momentum r₂ x V₂
-So, the component d₂ r₂  is not useful and you can delete  lines 417 to 420 , 405 to 408 , 392 to 400 , and lines 386 , 380 , 369 and 359
-Thus, 34 bytes may be saved.

-Lines 623 and 539 thru 561 are only useful if the eccentricity is exactly equal to 1  ( parabola )
-Practically, that seems highly improbable, so these lines may be deleted without taking many risks.
 
 

 
    ( 938 bytes / SIZE 042 )
 
 

 
Example:    Comet P2007/T2 Kowalski, on 2007/07/01 0h , 2007/07/05 0h , 2007/07/09 0h its position was:

                   2007/07/01 0h                   2007/07/05 0h                 2007/07/09 0h

      t                 2738                                   2742                                2746                               days
    RA         14h27m25s57                     14h16m33s94                   14h06m37s72
   Decl       -39°30'56"18                      -38°44'07"32                    -37°52'59"14
     X           -0.155835                           -0.222301                        -0.287788                   Astronomical Units
     Y            1.004614                             0.992089                         0.975105                     ------------------

-Store these 15 numbers into

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

  and  XEQ "GAUSS"   the successive  D₂-values are displayed ( the last displayed value is 0.593970 ) and eventually:

                      T      =  2817.895386               = R16                   --- execution time = 59 seconds ---
                      q      =   0.697257                    = R17
                      e      =   0.775182                    = R18
                       i      =   9°8769                       = R19
                 omega   = -1°4196 = 358°5804    = R20
               OMEGA =   3°8809                       = R21
                      n      =   0.180451                    = R22

-If we reduce these elements to the standard epoch J2000 ( cf for instance "Orbital Elements from one Equinox to another one" ) it yields:

                      i₂₀₀₀        =   9°8759
                 omega₂₀₀₀   = 358°5796
               OMEGA₂₀₀₀ =   3°7770

-In fact, the exact mean elements are

                           T       = 2818.01589 = 2007/09/19.01589
                           q        =  0.695805
                           e        =  0.774729
                           i         =  9°8974
                      omega    = 358°5346                     referred to J2000.0
                    OMEGA  =  4°0019
                          n         =  0.181564 deg/day

-The error is about 3 hours for the time of passage in perihelion but some of these results are relatively disappointing, especially OMEGA!
-However, we have found the elements of the osculating orbit which may be significantly different from the mean orbit.
-Like "OLBERS" , "GAUSS"  neglects the aberration, parallax and planetary perturbations.

Notes:

-The equation that is solved by lines 255 to 297 may have several solutions ( up to 3 ) and another evaluation a few days later may be useful...
-The formulas involve several cross products and dot products, therefore many bytes can be saved
  if you have CROSS and DOT micro-code routines ( cf for instance "A few M-Code Routines for the HP-41" ).
-Gauss' method is very sensitive to observational errors, so take the average of successive results to increase the accuracy.
-The program uses the ecliptic and equator of the date(s) and that could be criticized because this is not an inertial frame.
  However, the resulting errors are usually negligible - provided the intervals of time remain of the order of a few days.
-Of course, a better program should use a standard equinox - for instance J2000 - but in this case, the rectangular coordinates of the Sun
  will not verify Z ~ 0 and more data registers and more bytes would be required ( cf next paragraph ).

Improvements:
 

1°) Slightly more accurate results may be obtained if you

-replace lines 298 to 305 by

              RCL 38        /         RCL 28       ST* 34       RCL 35
              RCL 42       1           /                 ST* 40         +
              RCL 28       +        ST* 33         *

-replace line 282  by

             RCL 42         RCL 39
             RCL 28         ST* Y
                /                     +

-and add

             RCL 31          RCL 11        X^2            12               /                           after line 216
             RCL 32          RCL 01        ST* Y        RCL 41      STO 42
                  *                     -                   +                *

-These lines use Herrick-Gibbs formulae:

    c_i = A_i + ( B_i / r₂³ )( 1 + B₂ / r₂³ )  with  B₂ = ( k²/12 ) ( t₃ - t₁ )² ( 1 + A₁ A₃ )                             ( i = 1 , 3 )
 

2°) If you want to take the effects of aberration and light-time into account:

  add   FS?C 01        after line 321 ( STO 27 )
           GTO 16

  replace line 217 by  RCL 42

  and  add

        SF 01                LBL 16           /                      XEQ "SUN2"       RCL 40         XEQ "SUN2"            RCL 40        XEQ "SUN2"    after line 02
        CLX                  RCL 06          -                      P-R                       /                    P-R                            /                   P-R
        STO 00             RCL 00          STO 06           STO 09                -                    STO 04                      -                  STO 14
        STO 26             X#0?              365250           X<>Y                   STO 01         X<>Y                        STO 11        X<>Y
        STO 27             STO 42          STO 41           STO 10                RCL 41         STO 05                     RCL 41        STO 15
        SIGN                173.142          /                      RCL 01                 /                    RCL 11                     /
        STO 42             STO 40          STO 00           RCL 26                STO 00         RCL 27                     STO 00

 (  "SUN3" should produce more accurate results  - cf "Astronomical Ephemeris for the HP-41" §9 )

-Registers  R04-05-09-10-14-15  do not have to be initialized anymore.
 

     b)  Program#2
 

-If you want to use more accurate positions of the Sun and/or if you want to use the standard ecliptic J2000, you'll need more registers to store the Z_i
-"GAUSS+" uses the Herrick-Gibbs formulae mentioned above.
-Several M-code functions ( ST<>A ,  CROSS , DOT , ... ) are also used ( cf for instance "A few M-code Routines for the HP-41" ),
  but alternatives are suggested below.
 

Data Registers:              R00 & R19 to R49: 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   ( 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 ) - virtual for a hyperbolic orbit.

   >>>>  You have also:

                 R41 = Distance Earth-Comet at the 1st instant     R44 = Distance Sun-Comet at the 1st instant       = r₁
                 R00 = Distance Earth-Comet at the 2nd instant    R45 = (Distance Sun-Comet at the 2nd instant)^3 = r₂^3
                 R43 = Distance Earth-Comet at the 3rd instant     R46 = Distance Sun-Comet at the 3rd instant      = r₃

      The coordinates of the vector r₁ are stored in R31-R32-R33
      The coordinates of the vector r₂ are stored in R34-R35-R36
      The coordinates of the vector r₃ are stored in R37-R38-R39

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

-If you don't want to use M-code functions, add

        LBL 09      X<>Y        +               *                  after line 235  and replace  DOT  by  XEQ 09
        RCL 50     RCL 51     X<>Y        +
        *                *               RCL 52     RTN

  replace STO M  STO N  STO O  by  STO 50  STO 51  STO 52
    and     RCL M  RCL N  RCL O  by  RCL 50  RCL 51  RCL 52

-Use the "CROSS" routine listed in "Dot & Cross product for the HP-41" and follow the other suggestions below

-Lines 55 to 62 may be replaced by    34   RCL 39  RCL 38  RCL 37  XEQ "CROSS"
-Lines 68 to 71 may be replaced by    31   RCL 39  RCL 38  RCL 37  XEQ "CROSS"
-Lines 77 to 84 may be replaced by    31   RCL 36  RCL 35  RCL 34  XEQ "CROSS"

-Line 241 , this loop solves an equation by the secant method
-If the process seems to diverge or if it converges to a negative value:
-Stop the program, set flag F00 , store another guess into R00 and press  XEQ 10.

-Lines 300-301 may be replaced by   1   +   RCL 45
-Lines 376-384-501  may be replaced by  3  Y^X
-Lines 409 to 412 may be replaced by     34   RCL 52  RCL 51  RCL 50  XEQ 10
-Lines 455-480-509 are equivalent to  2   /
-Lines  534 to 537 may be replaced with   1   RCL Y   -   /   ST+ X   LN1+X   ENTER^   E^X-1   LASTX   CHS   E^X-1   -   2   /
 
 

 
     ( 872 bytes / SIZE 050 )
 
 

 
Example1:    Comet P2007/T2 Kowalski, on 2007/07/01 0h , 2007/07/05 0h , 2007/07/09 0h , astrometric coordinates ( /J2000 ):

                   2007/07/01 0h                   2007/07/05 0h                 2007/07/09 0h

      t                 2738                                   2742                                2746                               days
    RA       14h26m56s630                   14h16m05s582                 14h06m09s943
   Decl       -39°28'38"88                      -38°41'45"79                    -37°50'34"44
     X        -0.154038961                     -0.220524792                  -0.286041210                  Astronomical Units
     Y         1.004896850                       0.992492986                   0.975629006                   Astronomical Units
     Z        -0.000017928                     -0.000015437                  -0.000012870                   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

  SF 01  XEQ "GAUSS+"   the successive  D₂-values are displayed ( the last displayed value is 0.594663 ) and eventually:

                      T      =  2817.969991               = R19                   --- execution time = 65 seconds ---
                      q      =   0.696446                    = R20
                      e      =   0.774784                    = R21
                       i      =   9°8875                       = R22
                 omega   = -1°4533 = 358°5467    = R23
               OMEGA =   3°9160                       = R24
                      n      =   0.181247                    = R25

-The exact mean elements are:

                                     T       = 2818.01589 = 2007/09/19.01589
                                     q        =  0.695805
                                     e        =  0.774729
                                     i         =  9°8974
                                 omega    = 358°5346                     referred to J2000.0
                               OMEGA  =  4°0019
                                    n         =  0.181564 deg/day

Example2:     Comet C/2007 K3 Siding Spring on  2008/06/01  0h ,  2008/06/04   0h ,   2008/06/07   0h  TT , mean coordinates of the date:

        t             3074                                     3077                                            3080                                 days
      RA     22h03m09s301                   22h06m25s468                             22h09m28s952
     Decl       2°17'49"44                         3°11'59"67                                    4°05'25"16
       X        0.332127259                      0.283777164                                0.234699580                Astronomical Units
       Y        0.958175038                      0.974055899                                0.987437299                -------------------
       Z        0.000003049                       0.000002941                                0.000001543                -------------------

-Store these 18 numbers into

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

  CF 01  XEQ "GAUSS+"   the successive  D₂-values are displayed but ... the algorithm seems to converge to a negative value!
                                           So stop the program, set flag F00  ( SF 00 ), store another starting-value into R00 ( say,  2  STO 00 ) and  XEQ 10
                                           Now the process converges to 1.709255 and finally:
 

                      T      =  3033.66930                 = R19
                      q      =   2.050725                    = R20
                      e      =   1.001541                    = R21                         ( hyperbola )
                       i      =    16°2979                     = R22
                 omega   =    23°5733                     = R23
               OMEGA = -96°6300 = 263°3700   = R24

-After reducing these elements to the standard epoch J2000, it yields:

                      i₂₀₀₀       =  16°2979
                 omega₂₀₀₀   =  23°5772
               OMEGA₂₀₀₀ = 263°2485

-The exact mean elements are:

                   T       = 3033.66811 = 2008/04/21.66811
                   q        =  2.050848
                   e        =  1.001369
                   i         =  16°2998
               omega    =  23°5791           referred to J2000.0
             OMEGA  = 263°2551
 

Notes:

-You can add lines to take into account the case e = 1, but this doesn't seem very useful.
-Other lines may also be added to "GAUSS+" to take the effect of aberration and light-time into account:
-Use the distances Comet-Earth at the 3 instants in registers  R41  R00  R43.
 

References:

[1]  S. Bouiges - "Calcul Astronomique pour Amateurs" - Masson - ISBN  2-225-78265-2  ( in French )
[2]  Jean Meeus - "Astronomical Algorithms" - Willmann-Bell  -  ISBN 0-943396-35-2
[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