Comets

# 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  t1  the right ascension of the comet =  RA1  its declination = d1 and the rectangular ecliptic coordinates of the Sun are  X1 & Y1
-At an instant  t2  the right ascension of the comet =  RA2  its declination = d2 and the rectangular ecliptic coordinates of the Sun are  X2 & Y2
-At an instant  t3  the right ascension of the comet =  RA3  its declination = d3 and the rectangular ecliptic coordinates of the Sun are  X3 & Y3

-We must have   t1 <  t2 <  t3

>>> The results will be more accurate if  t3 - t2 is equal ( or nearly equal ) to   t2 - t1  ( this value should be of the order of a few days )

-Let

li  =  cos RAi cos di
mi = sin obl sin di + cos obl cos di sin RAi    where   obl is the obliquity of the ecliptic,   i = 1 , 2 , 3
ni = cos obl sin di - sin obl cos di sin RAi

-The ecliptic longitudes and latitudes of the comet are  Li = Atan mi/li  ( use R-P to get the correct quadrant ) and  bi = Asin ni
-Let Di = distance Earth-Comet at the instant ti , the Di 's are unknown but   Q = D3/D1 may be calculated by the approximate formula:

Q = ((t3 - t2)/(t2 - t1)) ( sin b1 / sin b3 ) { [ sin ( L2 - L'2 ) / tan b2 - sin ( L1 - L'2 ) / tan b1 ] / [ sin ( L3 - L'2 ) / tan b3 - sin ( L2 - L'2 ) / tan b2 ] }

where  L'2  is the ecliptic longitude of the Sun at the instant  t2

-"OLBERS" uses the equivalent formula:   Q = - ((t3 - t2)/(t2 - t1)) [ u1. ( u2 x R2 ) ] / [ u3. ( u2 x R2 ) ]       where  ui ( li , mi , ni )  and  R2 ( X2 , Y2 , 0 )

-The rectangular heliocentric coordinates of the comet at the instant   t1 &  t3  are then:

x1 = l1 D1 - X1               x3 = Q l3 D1 - X3
y1 = m1 D1 - Y1             y3 = Q m3 D1 - Y3
z1 = n1 D1                      z3 = Q n3 D1

-So the square of the distance  c2 = (x3 - x1)2 + (y3 - y1)2 + (z3 - z1)2    may be written under the form:   c2 = A D12 + B D1 + C

-On the other hand:

r12 = D12 - 2 D1 ( l1 X1 + m1 Y1 ) + X12 + Y12                             where  ri = radius vector of the comet at the instant  ti
r32 = Q2 D12 - 2.Q D1 ( l3 X3 + m3 Y3 ) + X32 + Y32

-We have also   c = ( r1 + r3 ) sin 2µ

with  µ = Asin { sqrt(2) Sin [ (1/3) Asin [ 3k ( t3 - t1 ) ( r1 + r3 ) -3/2 /sqrt(2) ] ] }     where  k = 0.01720209895 = Gaussian gravitational constant

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

-When D1 is found, we compute the quantities

s1 = x1 y3 - x3 y1
s2 = y1 z3 - y3 z1          and    S = ( s12 + s22 + s32 )1/2
s3 = x1 z3 - x3 z1

-The longitude of the ascending node is then   OMEGA = Atan2 s2/s3  ( use R-P to get the proper quadrant )
and the inclination                               i       = Acos s1/S

-    ( vi + omega ) may be obtained by

ri Sin ( vi + omega ) = zi / Sin i
ri Cos ( vi + omega ) = xi Cos OMEGA + yi Sin OMEGA

-The relation  v3 - v1 = Asin S/(r1 r3)  is also used - we assume that  v3 - v1   does not exceed  90°

-The perihelion distance  q = [ r1 r3 Sin2 (v3 - v1)/2 ] / [ r1 + r3 - 2(r1 r3)1/2 Cos (v3 - v1)/2 ]

-  (v3 + v1) is calculated with the R-P function from the relations

2.S Sin (v3 + v1)/2 = q ( r3 - r1 ) Cos (v3 - v1)/2
2.S Cos (v3 + v1)/2 = [ q ( r3 + r1 ) - r1 r3 ] Sin (v3 - v1)/2

whence  v1  whence  omega since we already know ( vi + omega )

-Finally    T = t1 - sqrt(2)/(3k) q3/2 ( 3.s + s3 )    with  s = Tan v1/2

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

•  R01 = t1               •  R06 = t2               •  R11 = t3                  expressed in days since 2000/01/01  0h TT
•  R02 = RA1           •  R07 = RA2           •  R12 = RA3             in hh.mnss
•  R03 = d1              •  R08 = d2              •  R13 = d3                 in °. ' ''
•  R04 = X1             •  R09 = X2             •  R14 = X3                 in Astronomical Units
•  R05 = Y1             •  R10 = Y2             •  R15 = Y3                 ---------------------

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

 01  LBL "OLBERS"   02  DEG   03  RCL 03   04  RCL 02   05  RCL 01   06  XEQ 01   07  STO 18   08  RDN   09  STO 17   10  X<>Y   11  STO 16   12  RCL 13   13  RCL 12   14  RCL 11   15  XEQ 01   16  STO 21   17  RDN   18  STO 20   19  X<>Y   20  STO 19   21  RCL 08   22  RCL 07   23  RCL 06   24  XEQ 01   25  STO 00   26  RCL 09   27  ST* Z   28  *   29  STO 22   30  X<> Z   31  RCL 10   32  ST* 00   33  *   34  -   35  STO 23   36  RCL 18   37  *   38  RCL 17   39  RCL 22   40  *   41  -   42  RCL 00   43  RCL 16   44  *   45  +   46  RCL 20   47  RCL 22   48  *   49  RCL 21   50  RCL 23   51  *   52  -   53  RCL 00   54  RCL 19   55  *   56  -   57  /   58  RCL 11   59  RCL 06   60  ST- Y   61  RCL 01   62  -   63  /   64  * 65  STO 24             66  1   67  STO 00   68  SF 00   69  GTO 16   70  LBL 01   71  2807410   72  /   73  23.439279   74  -   75  STO 00   76  RDN   77  HR   78  15   79  *   80  X<>Y   81  HR   82  STO 22   83  COS   84  P-R   85  X<>Y   86  RCL 00   87  X<>Y   88  P-R   89  X<>Y   90  X<> 00   91  RCL 22   92  SIN   93  P-R   94  ST+ 00   95  RDN   96  -   97  RCL 00   98  RTN   99  LBL 15  100  STO 28 101  .1 102  STO 27 103  ST+ 00  104  LBL 16 105  VIEW 00 106  RCL 00 107  RCL 04 108  RCL 16 109  * 110  RCL 05 111  RCL 17 112  * 113  + 114  ST+ X 115  - 116  * 117  RCL 04 118  X^2 119  + 120  RCL 05 121  X^2 122  + 123  SQRT 124  STO 22 125  CLX 126  RCL 24 127  * 128  RCL 14 129  RCL 19           130  * 131  RCL 15 132  RCL 20 133  * 134  + 135  ST+ X 136  - 137  * 138  RCL 14 139  X^2 140  + 141  RCL 15 142  X^2 143  + 144  SQRT 145  STO 23 146  RCL 22 147  + 148  STO 25 149  ENTER^ 150  SQRT 151  * 152  27.40389543 153  STO 29 154  * 155  RCL 11 156  RCL 01 157  - 158  X<>Y 159  / 160  ASIN 161  3 162  / 163  SIN 164  2 165  SQRT 166  * 167  ASIN 168  ST+ X 169  SIN 170  RCL 25 171  * 172  X^2 173  RCL 04 174  RCL 14 175  - 176  X^2 177  - 178  RCL 05 179  RCL 15 180  - 181  X^2 182  - 183  RCL 00 184  / 185  RCL 16 186  RCL 19 187  RCL 24 188  * 189  - 190  STO 25 191  RCL 04 192  RCL 14 193  - 194  * 195  ST+ X 196  + 197  RCL 17           198  RCL 20 199  RCL 24 200  * 201  - 202  STO 26 203  RCL 05 204  RCL 15 205  - 206  * 207  ST+ X 208  + 209  RCL 18 210  RCL 21 211  RCL 24 212  * 213  - 214  X^2 215  RCL 25 216  X^2 217  + 218  RCL 26 219  X^2 220  + 221  RCL 00 222  * 223  - 224  FS?C 00 225  GTO 15 226  ENTER^ 227  X<> 28 228  RCL Y 229  - 230  / 231  ST* 27 232  RCL 27 233  ST+ 00 234  ABS 235   E-6 236  XY 289  STO 20 290  RCL 24 291  R^ 292  / 293  ACOS 294  X<> 18 295  RCL 18 296  SIN 297  / 298  RCL 16 299  RCL 20 300  COS 301  * 302  RCL 17 303  RCL 20 304  SIN 305  * 306  + 307  R-P 308  X<>Y 309  X<> 19 310  RCL 22 311  RCL 23 312  * 313  STO 00 314  / 315  ASIN 316  2 317  / 318  STO 16 319  SIN 320  X^2 321  RCL 00           322  ST* Y 323  SQRT 324  RCL 16 325  COS 326  * 327  ST+ X 328  RCL 22 329  RCL 23 330  + 331  STO 24 332  X<>Y 333  - 334  / 335  STO 17 336  ST/ 00 337  RCL 16 338  X<>Y 339  P-R 340  RCL 23 341  RCL 22 342  - 343  * 344  X<>Y 345  RCL 24 346  RCL 00 347  - 348  * 349  R-P 350  CLX 351  RCL 16 352  - 353  ST- 19 354  2 355  / 356  TAN 357  ENTER^ 358  X^2 359  3 360  + 361  * 362  RCL 29 363  * 364  RCL 17 365  ST* Y 366  SQRT 367  * 368  RCL 01 369  X<>Y 370  - 371  STO 16 372  RCL 20 373  SIGN 374  RCL 19 375  RCL 18 376  RCL 17 377  RCL 16 378  CLD 379  END

( 562 bytes / SIZE 030 )

 STACK INPUTS OUTPUTS T / omega Z / i Y / q X / T L / OMEGA

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  D1-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:

i2000        = 117°6686
omega2000   = 233°6403
OMEGA2000 = 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:   s3 + 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+s2)

>>> 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 = lg = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
R08 = decl ( °. ' '' )                               R11 = bg = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
R09 = rT = 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

 01  LBL "PARM" 02  DEG 03  STO 00 04  RCL 01 05  - 06  54.80779086 07  / 08  RCL 02 09  ST/ Y 10  SQRT 11  / 12  SIGN 13  LASTX          14  ABS 15  ENTER^ 16  X^2 17  1 18  + 19  SQRT           20  + 21  3 22  1/X 23  Y^X 24  ENTER^ 25  1/X 26  - 27  * 28  ATAN 29  ST+ X 30  LASTX          31  X^2 32  RCL 02 33  ST* Y 34  + 35  XEQ "GEO"  36  END

( 62 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / Psi ( deg ) Z / rT ( AU ) Y / Decl ( °. ' '' ) X t ( days ) RA ( hh.mnss )

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         i2000      =   48°7131
q =  0.990662                                 omega2000   = 163°4788
e =  1                                            OMEGA2000 = 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     rT  =  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:

RA2000 =  16h20m08s
Decl2000 =  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 = lg = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
R08 = decl ( °. ' '' )                               R11 = bg = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
R09 = rT = 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 
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.

 01  LBL "ELLM"   02  DEG   03  STO 00   04  RCL 01   05  -   06  58.13244087   07  /   08  1   09  RCL 03   10  -   11  STO 09   12  RCL 02   13  /   14  ST* Y   15  SQRT   16  *   17  R-D   18  STO 07   19  LASTX   20  RCL 09   21  ST+ X   22  RCL 03   23  8   24  *   25  1 26  +   27  ST/ Z   28  /   29  STO 08           30  3   31  Y^X   32  RCL Y   33  X^2   34  +   35  SQRT   36  RCL Y   37  SIGN   38  *   39  +   40  SIGN   41  LASTX   42  ABS   43  3   44  1/X   45  Y^X   46  *   47  RCL 08   48  2   49  /   50  - 51  STO Y   52  5   53  Y^X   54  .078   55  *   56  RCL 03           57  1   58  +   59  /   60  -   61  3   62  RCL Y   63  ST+ X   64  X^2   65  -   66  *   67  RCL 03   68  *   69  R-D   70  RCL 07   71  +   72  STO 08   73  LBL 01   74  RCL 08   75  SIN 76  RCL 03   77  R-D   78  *   79  RCL 07   80  +   81  1   82  RCL 08   83  ST- Z   84  COS   85  RCL 03           86  *   87  -   88  /   89  ST+ 08   90  ABS   91   E-6   92  XY 108  ST+ X 109  ENTER^ 110  COS 111  RCL 03         112  ST* Y 113  1 114  ST+ Z 115  + 116  RCL 02 117  * 118  X<>Y 119  / 120  XEQ "GEO" 121  END

( 166 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / Psi ( deg ) Z / rT ( AU ) Y / Decl ( °. ' '' ) X t ( days ) RA ( hh.mnss )

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                 i2000      =  105°063204
q =  3.432968                                 omega2000   = 337°140230
e =  0.984585                               OMEGA2000 = 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     rT  =  4.4261 AU
RDN    Psi  =  38°52

-Referred to J2000 , it yields:

RA2000 =  16h20m08s
Decl2000 =  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 = lg = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
R08 = decl ( °. ' '' )                               R11 = bg = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
R09 = rT = 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.

 01  LBL "HYPM" 02  DEG 03  STO 00 04  RCL 01 05  - 06  58.13244087 07  / 08  RCL 03 09  1 10  - 11  STO 10 12  RCL 02 13  / 14  ST* Y 15  SQRT 16  * 17  STO 08 18  ABS 19  X#0? 20  LN 21  ABS 22  RCL 08         23  ABS 24  X>Y? 25  X<>Y 26  LASTX 27  SIGN 28  * 29  STO 07 30  LBL 01 31  RCL 07 32  ENTER^ 33  ST+ Y 34  E^X-1 35  STO 09 36  LASTX 37  CHS 38  E^X-1 39  ST+ 09 40  - 41  RCL 03         42  ST* 09 43  * 44  - 45  RCL 08 46  ST+ X 47  + 48  RCL 09 49  RCL 10 50  ST+ X 51  + 52  / 53  ST+ 07 54  ABS 55   E-6 56  XY 87  / 88  XEQ "GEO" 89  END

( 126 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / Psi ( deg ) Z / rT ( AU ) Y / Decl ( °. ' '' ) X t ( days ) RA ( hh.mnss )

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                i2000      = 117°649041
q =  0.969480                                          omega2000   = 233°671201
e =  1.000785                                        OMEGA2000 = 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     rT  =  1.5825 AU                                                                         RA2000 =  17h02m13s
RDN    Psi  =  39°16                                                                                Decl2000 =  -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" ,  rT  =  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 = lg = geocentric longitude ( in degrees )             R12 = r = radius vector ( AU )
R08 = decl ( °. ' '' )                               R11 = bg = geocentric latitude ( in degrees )               R13 = L = heliocentric longitude ( in degrees )
R09 = rT = 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.

 01  LBL "GEO" 02  STO 12        03  CLX 04  RCL 05 05  + 06  ENTER^ 07  SIN 08  RCL 04 09  SIN 10  * 11  ASIN 12  STO 14 13  X<>Y 14  1 15  P-R 16  X<>Y 17  RCL 04 18  COS 19  * 20  X<>Y 21  R-P 22  CLX 23  RCL 06 24  + 25  STO 13 26  365250 27  ST/ 00 28  XEQ "SUN2" 29  X<>Y 30  STO 07         31  X<>Y 32  P-R 33  STO 15 34  X<>Y 35  STO 16  36  RCL 14 37  RCL 13 38  RCL 12         39  XEQ "S-R" 40  RCL 16 41  ST+ Z 42  CLX 43  RCL 15 44  + 45  XEQ "R-S" 46  STO 09 47  CLX 48  .1301 49  RCL 00 50  * 51  23.43928 52  - 53  X<> Z 54  STO 11         55  X<>Y 56  ST- 07 57  STO 10 58  XEQ "EE" 59  15 60  / 61  24 62  MOD 63  HMS 64  RCL 07 65  COS 66  RCL 11         67  COS 68  * 69  ACOS 70  RCL 09 71  R^ 72  HMS 73  STO 08 74  R^ 75  STO 07 76  END

( 124 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / Psi ( deg ) Z / rT ( AU ) Y v ( deg ) Decl ( °. ' '' ) X r ( AU ) RA ( hh.mnss )

r = radius vector       RA = right-ascension     rT = 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  t1  the right ascension of the comet =  RA1  its declination = d1 and the rectangular ecliptic coordinates of the Sun are  X1 & Y1
-At an instant  t2  the right ascension of the comet =  RA2  its declination = d2 and the rectangular ecliptic coordinates of the Sun are  X2 & Y2
-At an instant  t3  the right ascension of the comet =  RA3  its declination = d3 and the rectangular ecliptic coordinates of the Sun are  X3 & Y3

-We must have   t1 <  t2 <  t3

>>> The results will be more accurate if  t3 - t2 is equal ( or nearly equal ) to   t2 - t1  ( this value should be of the order of a few days )

-Let

li  =  cos RAi cos di
mi = sin obl sin di + cos obl cos di sin RAi    where   obl is the obliquity of the ecliptic,   i = 1 , 2 , 3
ni = cos obl sin di - sin obl cos di sin RAi

-Let Di = distance Earth-Comet at the instant ti , ui ( li , mi , ni ) , -Ri = position vector of the Earth , ri = position vector of the Comet

ri = Di ui  - Ri        (  || ui || = 1  )

-If we neglect the planetary perturbations, the 3 position vectors ri  lie in the same plane, so there are coefficients c1 , c3 such that

r2 = c1 r1 + c3 r3   and likewise for the velocity at the second instant:        V2 = -d1 r1 + d2 r2 + d3 r3

-One can prove that, approximately:     ci = Ai + Bi / r23   ( i = 1 , 3 )  and   di = Gi + Hi / ri3  ( i = 1 , 2 , 3 )

where    A1 = (t3 - t2)/(t3 - t1)  ,   B1 = k2A1 [ (t3 - t1)2 - (t3 - t2)2 ]/6      with   k =  0.01720209895 = Gaussian gravitational constant
A3 = (t2 - t1)/(t3 - t1)  ,   B3 = k2A3 [ (t3 - t1)2 - (t2 - t1)2 ]/6

and      G1 = (t3 - t2)/[(t3 - t1)(t2 - t1)]       G3 = (t2 - t1)/[(t3 - t1)(t3 - t2)]        G2 = G1 - G3
H1 = k2(t3 - t2)/12                        H3 = k2(t2 - t1)/12                         H2 = H1 - H3

-Let   A( A1 X1 - X2 + A3 X3 , A1 Y1 - Y2 + A3 Y3 ,  A1 Z1 - Z2 + A3 Z3 )
and   B( B1 X1 + B3 X3 , B1 Y1 + B3 Y3 ,  B1 Z1 + B3 Z3 )                               ( all the Zi's ~ 0 if we use the rectangular ecliptic coordinates of the Sun )

P1 = A.( u2 x u3 )/D          Q1 = B.( u2 x u3 )/D
P2 = A.( u1 x u3 )/D          Q2 = B.( u1 x u3 )/D          where      D = u1.( u2 x u3 )           ( . = dot product , x = cross product )
P3 = A.( u1 x u2 )/D          Q3 = B.( u1 x u2 )/D

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

D2 = P2 + Q2 / r23   with    r2 =  D22 - 2 D2 ( l2 X2 + m2 Y2 ) + X22 + Y22

-Then,  D1 = ( P1 + Q1 / r23 )/( A1 + B1 / r23 )  and  D3 = ( P3 + Q3 / r23 )/( A3 + B3 / r23 )

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

xi = li Di - Xi
yi = mi Di - Yi        i = 1 , 2 , 3         so we know the 3 vectors ri's   and the velocity V2 = -d1 r1 + d2r2 + d3 r3   may be computed
zi = ni Di

-We calculate the cross-product  h = r2 x V2 and then:

-The parameter  p = h2/k2
-The inclination  i = Acos hz/h = Atan2 (hx2+hy2)1/2/hz                          you can also use C = r1 x r3  instead of  h  to get i and OMEGA

-The longitude of the ascending node  OMEGA = Atan2 hx/(-hy)

-    ( vi + omega ) is obtained by

ri Sin ( vi + omega ) = zi / Sin i
ri Cos ( vi + omega ) = xi Cos OMEGA + yi Sin OMEGA        whence   (v3 - v1)/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 (v3 + v1)/2 = ( p/r1 - p/r3 ) / Sin (v3 - v1)/2                   give  e and (v3 + v1)/2  whence v1 and omega
2e Cos (v3 + v1)/2 = ( p/r1 + p/r3 - 2 ) / Cos (v3 - v1)/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 = t1               •  R06 = t2               •  R11 = t3                  expressed in days since 2000/01/01  0h TT
•  R02 = RA1           •  R07 = RA2           •  R12 = RA3             in hh.mnss
•  R03 = d1              •  R08 = d2              •  R13 = d3                 in °. ' ''
•  R04 = X1             •  R09 = X2             •  R14 = X3                 in Astronomical Units
•  R05 = Y1             •  R10 = Y2             •  R15 = Y3                 ---------------------

>>>> 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 r2 x V2
-So, the component d2 r2  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.

 01  LBL "GAUSS"   02  DEG   03  RCL 03   04  RCL 02   05  RCL 01   06  XEQ 01   07  STO 18   08  RDN   09  STO 17   10  X<>Y   11  STO 16   12  RCL 08   13  RCL 07   14  RCL 06   15  XEQ 01   16  STO 21   17  RDN   18  STO 20   19  X<>Y   20  STO 19   21  RCL 13   22  RCL 12   23  RCL 11   24  XEQ 01   25  STO 24   26  RDN   27  STO 23   28  RCL 19   29  *   30  X<>Y   31  STO 22   32  RCL 20   33  *   34  -   35  STO 00   36  RCL 21   37  RCL 22   38  *   39  RCL 19   40  RCL 24   41  *   42  -   43  STO 26   44  RCL 20   45  RCL 24   46  *   47  RCL 21   48  RCL 23   49  *   50  -   51  STO 25   52  RCL 17   53  RCL 24   54  *   55  RCL 18   56  RCL 23   57  *   58  -   59  STO 27   60  RCL 18   61  RCL 22   62  *   63  RCL 16   64  RCL 24   65  *   66  -   67  STO 28   68  RCL 17   69  RCL 21   70  *   71  RCL 18   72  RCL 20   73  *   74  -   75  STO 29   76  RCL 18   77  RCL 19   78  *   79  RCL 16   80  RCL 21   81  *   82  -   83  STO 30   84  RCL 11   85  RCL 06   86  -   87  STO 31   88  X^2   89  RCL 06   90  RCL 01   91  -   92  STO 32   93  X^2   94  RCL 11   95  RCL 01   96  -   97  ST/ 31   98  ST/ 32   99  X^2 100  ST- Z 101  X<>Y 102  - 103  RCL 32 104  * 105  STO 34 106  X<>Y 107  CHS 108  RCL 31 109  * 110  3379.380681 111  STO 41 112  6 113  * 114  ST/ 34 115  / 116  STO 33 117  RCL 04 118  * 119  RCL 14 120  RCL 34 121  * 122  + 123  STO 39 124  RCL 25 125  * 126  RCL 05 127  RCL 33 128  * 129  RCL 15 130  RCL 34 131  * 132  + 133  STO 40 134  RCL 26 135  * 136  + 137  RCL 16 138  RCL 25 139  * 140  RCL 17 141  RCL 26 142  * 143  + 144  RCL 18 145  RCL 00 146  * 147  + 148  STO 00 149  / 150  STO 38 151  RCL 27 152  RCL 39 153  * 154  RCL 28 155  RCL 40 156  * 157  + 158  RCL 00 159  / 160  X<> 39 161  RCL 29 162  * 163  RCL 30 164  RCL 40 165  * 166  + 167  RCL 00 168  / 169  STO 40 170  RCL 04 171  RCL 31 172  * 173  RCL 14 174  RCL 32 175  * 176  + 177  RCL 09 178  - 179  STO 36 180  RCL 25 181  * 182  RCL 05 183  RCL 31 184  * 185  RCL 15 186  RCL 32 187  * 188  + 189  RCL 10 190  - 191  STO 37 192  RCL 26 193  * 194  + 195  RCL 00 196  / 197  STO 35 198  RCL 27 199  RCL 36 200  * 201  RCL 28 202  RCL 37 203  * 204  + 205  RCL 00 206  / 207  X<> 36 208  RCL 29 209  * 210  RCL 30 211  RCL 37 212  * 213  + 214  RCL 00          215  / 216  STO 37 217  1 218  STO 00 219  SF 00 220  GTO 10 221  LBL 01 222  2807410 223  / 224  23.439279 225  - 226  STO 00 227  RDN 228  HR 229  15 230  * 231  X<>Y 232  HR 233  STO 25 234  COS 235  P-R 236  X<>Y 237  RCL 00 238  X<>Y 239  P-R 240  X<>Y 241  X<> 00 242  RCL 25 243  SIN 244  P-R 245  ST+ 00 246  RDN 247  - 248  RCL 00 249  RTN 250  LBL 09 251  STO 26 252  .1 253  STO 27 254  ST+ 00 255  LBL 10 256  VIEW 00 257  RCL 00 258  RCL 09 259  RCL 19 260  * 261  RCL 10 262  RCL 20 263  * 264  + 265  ST+ X 266  - 267  * 268  RCL 09 269  X^2 270  + 271  RCL 10 272  X^2 273  + 274  ENTER^ 275  SQRT 276  * 277  STO 28 278  X<>Y 279  RCL 36 280  - 281  * 282  RCL 39 283  - 284  FS?C 00 285  GTO 09 286  ENTER^ 287  X<> 26 288  RCL Y 289  - 290  / 291  ST* 27 292  RCL 27 293  ST+ 00 294  ABS 295   E-7 296  XY 461  STO 21 462  CLX 463  RCL 39 464  R-P 465  X<>Y 466  STO 19 467  RCL 18 468  X<>Y 469  SIN 470  / 471  RCL 16 472  RCL 21 473  COS 474  * 475  RCL 17 476  RCL 21 477  SIN 478  * 479  + 480  R-P 481  X<>Y 482  STO 20 483  RCL 24 484  RCL 19 485  SIN 486  / 487  RCL 22 488  RCL 21 489  COS 490  * 491  RCL 23 492  RCL 21 493  SIN 494  * 495  + 496  R-P 497  CLX 498  RCL 20 499  - 500  2 501  / 502  STO 22 503  RCL 40 504  RCL 29 505  / 506  STO Y 507  RCL 40 508  RCL 30 509  / 510  ST+ Z 511  - 512  RCL 22 513  SIN 514  / 515  X<>Y 516  2 517  - 518  RCL 22 519  COS 520  / 521  R-P 522  2 523  / 524  STO 18 525  CLX 526  RCL 22 527  - 528  ST- 20 529  STO 22 530  RCL 40          531  1 532  RCL 18 533  + 534  / 535  STO 17 536  1 537  RCL 18 538  - 539  X#0? 540  GTO 00 541  X<> 22 542  2 543  / 544  TAN 545  ENTER^ 546  X^2 547  3 548  + 549  * 550  RCL 17 551  3 552  Y^X 553  ST+ X 554  9 555  / 556  RCL 41 557  * 558  SQRT 559  * 560  GTO 04 561  LBL 00 562  X<0? 563  SF 00 564  ABS 565  STO 33 566  / 567  3 568  Y^X 569  RCL 41 570  * 571  SQRT 572  1/X 573  FC? 00 574  R-D 575  X<> 22 576  2 577  / 578  RCL 33 579  1 580  RCL 18 581  + 582  / 583  SQRT 584  FS?C 00 585  GTO 02 586  P-R 587  LASTX 588  / 589  R-P 590  X<>Y 591  ST+ X 592  ENTER^ 593  SIN 594  RCL 18 595  R-D 596  * 597  GTO 03 598  LBL 02 599  X<>Y 600  TAN 601  * 602  1 603  RCL Y 604  - 605  / 606  ST+ X 607  LN1+X 608  ENTER^ 609  E^X-1 610  LASTX 611  CHS 612  E^X-1 613  - 614  2 615  / 616  RCL 18 617  * 618  X<>Y 619  LBL 03 620  - 621  RCL 22 622  / 623  LBL 04 624  RCL 01 625  X<>Y 626  - 627  STO 16 628  CLD 629  END

( 938 bytes / SIZE 042 )

 STACK INPUT OUTPUT X / T

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  D2-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:

i2000        =   9°8759
omega2000   = 358°5796
OMEGA2000 =   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:

ci = Ai + ( Bi / r23 )( 1 + B2 / r23 )  with  B2 = ( k2/12 ) ( t3 - t1 )2 ( 1 + A1 A3 )                             ( 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 Zi
-"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 = t1               •  R07 = t2               •  R13 = t3                  expressed in days since 2000/01/01  0h TT
•  R02 = RA1           •  R08 = RA2          •  R14 = RA3              in hh.mnss
•  R03 = d1              •  R09 = d2              •  R15 = d3                 in °. ' ''
•  R04 = X1             •  R10 = X2             •  R16 = X3                 in Astronomical Units          rectangular ecliptic
•  R05 = Y1             •  R11 = Y2             •  R17 = Y3                 ---------------------              coordinates
•  R06 = Z1             •  R12 = Z2              •  R18 = Z3                 ---------------------               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       = r1
R00 = Distance Earth-Comet at the 2nd instant    R45 = (Distance Sun-Comet at the 2nd instant)^3 = r2^3
R43 = Distance Earth-Comet at the 3rd instant     R46 = Distance Sun-Comet at the 3rd instant      = r3

The coordinates of the vector r1 are stored in R31-R32-R33
The coordinates of the vector r2 are stored in R34-R35-R36
The coordinates of the vector r3 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   /

 01 LBL "GAUSS+"   02  DEG   03  15   04  STO 19   05  39   06  STO 20   07  LBL 01   08  RCL IND 19   09  DSE 19   10  RCL IND 19   11  DSE 19   12  FS? 01   13  0   14  FC? 01   15  RCL IND 19   16  2807410   17  /   18  23.439279   19  -   20  STO 00   21  RDN   22  HR   23  15   24  *   25  X<>Y   26  HR   27  STO 21   28  COS   29  P-R   30  X<>Y   31  RCL 00   32  X<>Y   33  P-R   34  X<>Y   35  X<> 00   36  RCL 21   37  SIN   38  P-R   39  ST+ 00   40  RDN   41  -   42  RCL 00   43  STO IND 20   44  DSE 20   45  RDN   46  STO IND 20   47  DSE 20   48  X<>Y   49  STO IND 20   50  3   51  ST- 19   52  DSE 20   53  DSE 19   54  GTO 01   55  RCL 39   56  RCL 38   57  RCL 37   58  ST<>A   59  RCL 36   60  RCL 35   61  RCL 34   62  CROSS   63  STO 40   64  RDN   65  STO 41   66  X<>Y   67  STO 42   68  RCL 33   69  RCL 32   70  RCL 31   71  CROSS   72  STO 43   73  RDN   74  STO 44   75  X<>Y   76  STO 45   77  RCL 36   78  RCL 35   79  RCL 34   80  ST<>A   81  RCL 33   82  RCL 32   83  RCL 31   84  CROSS   85  STO 46   86  RCL 37   87  *   88  X<>Y   89  STO 47   90  RCL 38   91  *   92  + 93  X<>Y   94  STO 48   95  RCL 39             96  *   97  +   98  STO 00   99  RCL 13 100  RCL 07 101  - 102  STO 19 103  X^2 104  RCL 07 105  RCL 01 106  - 107  STO 20 108  X^2 109  RCL 13 110  RCL 01 111  - 112  ST/ 19 113  ST/ 20 114  X^2 115  ST- Z 116  X<>Y 117  - 118  RCL 20 119  * 120  STO 22 121  X<>Y 122  CHS 123  RCL 19 124  * 125  3379.380681 126  STO 30 127  6 128  * 129  ST/ 22 130  / 131  STO 21 132  RCL 04 133  * 134  RCL 16 135  RCL 22 136  * 137  + 138  STO M 139  RCL 05 140  RCL 21 141  * 142  RCL 17 143  RCL 22 144  * 145  + 146  STO N 147  RCL 06 148  RCL 21 149  * 150  RCL 18 151  RCL 22 152  * 153  + 154  RCL 00 155  ST/ M 156  ST/ N 157  / 158  STO O 159  29 160  STO 23 161  XEQ 02 162  RCL 04 163  RCL 19 164  * 165  RCL 10 166  - 167  RCL 16 168  RCL 20 169  * 170  + 171  STO M 172  RCL 05 173  RCL 19 174  * 175  RCL 11 176  - 177  RCL 17 178  RCL 20 179  * 180  + 181  STO N 182  RCL 06 183  RCL 19 184  * 185  RCL 12 186  - 187  RCL 18           188  RCL 20 189  * 190  + 191  RCL 00 192  ST/ M 193  ST/ N 194  / 195  STO O 196  DSE 23 197  XEQ 02 198  RCL 19 199  RCL 20 200  * 201  RCL 13 202  RCL 01 203  - 204  X^2 205  ST* Y 206  + 207  12 208  RCL 30 209  * 210  STO 49 211  / 212  STO 23 213  1 214  STO 00 215  SF 00 216  GTO 10 217  LBL 02 218  RCL 48 219  RCL 47 220  RCL 46 221  DOT 222  STO IND 23 223  DSE 23 224  RCL 45 225  RCL 44 226  RCL 43 227  DOT 228  STO IND 23 229  DSE 23 230  RCL 42 231  RCL 41 232  RCL 40 233  DOT 234  STO IND 23 235  RTN 236  LBL 03 237  STO 40 238  .1 239  STO 41 240  ST+ 00 241  LBL 10 242  VIEW 00 243  RCL 00 244  RCL 10 245  RCL 34 246  *  247  RCL 11 248  RCL 35 249  * 250  + 251  RCL 12 252  RCL 36 253  * 254  + 255  ST+ X 256  - 257  * 258  RCL 10 259  X^2 260  + 261  RCL 11 262  X^2 263  + 264  RCL 12 265  X^2 266  + 267  ENTER^ 268  SQRT 269  * 270  STO 45 271  X<>Y 272  RCL 25 273  - 274  * 275  RCL 23 276  RCL 45 277  / 278  RCL 28           279  ST* Y 280  + 281  - 282  FS?C 00 283  GTO 03 284  ENTER^ 285  X<> 40 286  RCL Y 287  - 288  / 289  ST* 41 290  RCL 41 291  ST+ 00 292  ABS 293   E-7 294  XY 414  CHS 415  R-P 416  RCL Z 417  R-P 418  STO 20 419  CLX 420  RCL 33 421  X<>Y 422  STO 22 423  SIN 424  / 425  X<>Y 426  STO 24 427  COS 428  RCL 31 429  * 430  RCL 32 431  RCL 24 432  SIN 433  * 434  + 435  R-P 436  X<>Y 437  STO 23 438  RCL 39 439  RCL 22 440  SIN 441  / 442  RCL 37 443  RCL 24 444  COS 445  * 446  RCL 38 447  RCL 24 448  SIN 449  * 450  + 451  R-P 452  CLX 453  RCL 23 454  - 455  X/2 456  STO 29 457  RCL 20 458  X^2 459  RCL 30 460  * 461  STO 20 462  RCL 44           463  / 464  STO Y 465  RCL 20 466  RCL 46 467  / 468  ST+ Z 469  - 470  RCL 29 471  SIN 472  / 473  X<>Y 474  2 475  - 476  RCL 29 477  COS 478  / 479  R-P 480  X/2 481  STO 21 482  CLX 483  RCL 29 484  - 485  ST- 23 486  STO 25 487  RCL 20 488  1 489  RCL 21 490  + 491  / 492  STO 20 493  1 494  RCL 21 495  - 496  X<0? 497  SF 00 498  ABS 499  STO 26 500  / 501  X^3  502  RCL 30 503  * 504  SQRT 505  1/X 506  FC? 00 507  R-D 508  X<> 25 509  X/2 510  RCL 26 511  1 512  RCL 21 513  + 514  / 515  SQRT 516  FS?C 00 517  GTO 04 518  P-R 519  LASTX 520  / 521  R-P 522  X<>Y 523  ST+ X 524  ENTER^ 525  SIN 526  RCL 21 527  R-D 528  * 529  GTO 05 530  LBL 04 531  X<>Y 532  TAN 533  * 534  ATANH 535  ST+ X 536  ENTER^ 537  SINH 538  RCL 21 539  *  540  X<>Y 541  LBL 05 542  - 543  RCL 25 544  / 545  RCL 01 546  X<>Y 547  - 548  STO 19 549  CLA 550  CLD 551  END

( 872 bytes / SIZE 050 )

 STACK INPUT OUTPUT X / T

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  D2-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  D2-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:

i2000       =  16°2979
omega2000   =  23°5772
OMEGA2000 = 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:

  S. Bouiges - "Calcul Astronomique pour Amateurs" - Masson - ISBN  2-225-78265-2  ( in French )
  Jean Meeus - "Astronomical Algorithms" - Willmann-Bell  -  ISBN 0-943396-35-2
  Robin M. Green - "Spherical Astronomy" - Cambridge University Press - ISBN  0-521-31779-7
  John D. Anderson - JPL Technical Report n° 32-497
  Dan Boulet - "Methods of Orbit Determination for the Micro Computer" - Willmann-Bell  -  ISBN 0-943396-34-4