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
 
 

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  X<Y? 237  GTO 16 238  RCL 00 239  ST* 16 240  ST* 17 241  ST* 18 242  RCL 24 243  * 244  ST* 19 245  ST* 20 246  ST* 21 247  RCL 04 248  ST- 16 249  RCL 05 250  ST- 17 251  RCL 14 252  ST- 19 253  RCL 15 254  ST- 20 255  RCL 16 256  RCL 21

257  * 258  RCL 18           259  RCL 19 260  * 261  - 262  STO 26 263  X^2 264  RCL 17 265  ST* 19 266  RCL 21 267  * 268  RCL 18 269  RCL 20 270  * 271  - 272  STO 25 273  X^2 274  + 275  RCL 16 276  RCL 20 277  * 278  RCL 19 279  - 280  STO 24 281  X^2 282  + 283  SQRT 284  STO 19 285  RCL 25 286  RCL 26 287  R-P 288  X<>Y 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  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
 
 

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

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  X<Y?   93  GTO 01   94  RCL 08   95  2   96  /   97  1   98  RCL 03   99  + 100  RCL 09

101  / 102  SQRT 103  P-R 104  LASTX 105  / 106  R-P 107  X<>Y 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                 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.
 
 

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  X<Y? 57  GTO 01 58  RCL 07         59  SIGN 60  LASTX 61  ABS 62  E^X-1 63  RCL X 64  2 65  + 66  / 67  * 68  1 69  RCL 03 70  + 71  RCL 10 72  /

73  SQRT 74  * 75  ATAN 76  ST+ X 77  ENTER^ 78  COS 79  RCL 03         80  ST* Y 81  1 82  ST+ Z 83  + 84  RCL 02 85  * 86  X<>Y 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                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.
 
 

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

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  X<Y? 297  GTO 10 298  RCL 35 299  RCL 38 300  RCL 28 301  ST/ 33 302  ST/ 34 303  ST/ 40 304  / 305  + 306  RCL 31 307  RCL 33 308  + 309  / 310  STO 26  311  ST* 16  312  ST* 17 313  ST* 18 314  RCL 37 315  RCL 40

316  + 317  RCL 32 318  RCL 34 319  + 320  / 321  STO 27          322  ST* 22 323  ST* 23 324  ST* 24 325  RCL 04 326  ST- 16 327  RCL 05 328  ST- 17 329  RCL 14 330  ST- 22 331  RCL 15 332  ST- 23 333  RCL 16 334  RCL 17 335  R-P 336  RCL 18 337  R-P 338  STO 29 339  RCL 22 340  RCL 23 341  R-P 342  RCL 24 343  R-P 344  STO 30 345  RCL 00 346  ST* 19 347  ST* 20 348  ST* 21 349  RCL 09 350  ST- 19 351  RCL 10 352  ST- 20 353  RCL 31 354  RCL 06 355  RCL 01 356  -  357  STO 35 358  /  359  STO 25 360  RCL 11 361  RCL 06 362  - 363  ST/ 32 364  12 365  RCL 41 366  * 367  ST/ 35 368  / 369  STO 34 370  RCL 29 371  3 372  Y^X 373  / 374  + 375  CHS 376  STO 37 377  RCL 16 378  * 379  RCL 35 380  ST- 34 381  RCL 30 382  3 383  Y^X 384  / 385  RCL 32 386  ST- 25 387  + 388  STO 39 389  RCL 22 390  * 391  + 392  RCL 25 393  RCL 34 394  RCL 28 395  / 396  + 397  STO 38 398  RCL 19 399  * 400  + 401  STO 34 402  RCL 17 403  RCL 37 404  * 405  RCL 20 406  RCL 38 407  * 408  + 409  RCL 23 410  RCL 39 411  * 412  + 413  STO 35 414  RCL 18 415  RCL 37 416  * 417  RCL 21 418  RCL 38 419  * 420  +

421  RCL 24 422  RCL 39 423  * 424  + 425  STO 36 426  RCL 20          427  * 428  RCL 21 429  RCL 35 430  * 431  - 432  STO 37 433  X^2 434  RCL 19 435  RCL 36 436  * 437  RCL 21 438  RCL 34 439  * 440  - 441  STO 38 442  X^2 443  + 444  RCL 19 445  RCL 35 446  * 447  RCL 20 448  RCL 34 449  * 450  - 451  STO 39 452  X^2 453  + 454  RCL 41 455  * 456  STO 40 457  RCL 37 458  RCL 38 459  R-P 460  X<>Y 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  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   /
 
 

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  X<Y? 295  GTO 10 296  RCL 27 297  RCL 23 298  RCL 45 299  / 300  X+1 301  LASTX 302  / 303  ST* 21 304  ST* 22 305  ST* 29 306  * 307  RCL 24 308  + 309  RCL 19 310  RCL 21 311  + 312  / 313  STO 41 314  ST* 31 315  ST* 32 316  ST* 33 317  RCL 26 318  RCL 29 319  + 320  RCL 20 321  RCL 22 322  + 323  / 324  STO 43 325  ST* 37 326  ST* 38 327  ST* 39 328  RCL 04 329  ST- 31 330  RCL 05 331  ST- 32 332  RCL 06 333  ST- 33 334  RCL 16 335  ST- 37 336  RCL 17 337  ST- 38 338  RCL 18 339  ST- 39 340  RCL 31 341  RCL 32 342  R-P 343  RCL 33 344  R-P 345  STO 44 346  RCL 37 347  RCL 38 348  R-P 349  RCL 39 350  R-P 351  STO 46 352  RCL 00 353  ST* 34 354  ST* 35 355  ST* 36 356  RCL 10 357  ST- 34 358  RCL 11 359  ST- 35 360  RCL 12 361  ST- 36 362  RCL 20 363  RCL 13 364  RCL 07 365  - 366  STO 47 367  / 368  RCL 07

369  RCL 01           370  - 371  ST/ 19 372  RCL 49 373  ST/ 47 374  / 375  RCL 46 376  X^3 377  / 378  + 379  STO 48 380  RCL 37 381  * 382  RCL 47 383  RCL 44 384  X^3 385  / 386  RCL 19 387  + 388  STO 47 389  RCL 31 390  * 391  - 392  STO M 393  RCL 38 394  RCL 48 395  * 396  RCL 32 397  RCL 47 398  * 399  - 400  STO N 401  RCL 39 402  RCL 48 403  * 404  RCL 33 405  RCL 47 406  * 407  - 408  STO O 409  RCL 36 410  RCL 35 411  RCL 34 412  CROSS 413  X<>Y 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  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