# hp41programs

Orbit Determination

# Orbit Determination for the HP-41

Overview

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

0°)  3 Subroutines

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

a) Rectangular-Spherical conversion:

x = r cos b cos l
y = r cos b sin l
z = r sin b

where    x , y , z = rectangular coordinates,    r  = ( x2 + y2 + z2 )1/2  ,   l  = longitude ,  b = latitude

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

 01  LBL "R-S"  02  R-P  03  X<>Y  04  RDN  05  R-P  06  R^  07  X<>Y  08  END

( 16 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T T T Z z b Y y l X x r L / (x2+y2)1/2

Example:     x = 3 ; y = 4 ; z = -7   find the spherical coordinates ( in DEG mode )

-7  ENTER^
4  ENTER^
3  XEQ "R-S"    r  =  8.602325267
RDN     = 53.13010235°
RDN     b = -54.46232221°

b) Spherical-Rectangular conversion:

 01  LBL "S-R"  02  X<>Y  03  RDN  04  P-R  05  R^  06  X<>Y  07  P-R  08  END

( 16 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T T T Z b z Y l y X r x L / (x2+y2)1/2

Example:     r = 10 ; l = 124° ; b = 37°   find    x ; y ; z   ( in DEG mode )

37   ENTER^
124  ENTER^
10  XEQ "S-R"   x = -4.465913097
RDN   y =  6.620988446
RDN   z =  6.018150232

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

c) A very useful subroutine: "EE"

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

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

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

 01  LBL "EE"  02  1  03  XEQ "S-R"  04  RDN  05  R-P  06  X<> Z  07  ST- Y  08  X<> Z  09  P-R  10  R^  11  XEQ "R-S"  12  RDN  13  END

( 31 bytes / SIZE 000 )

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

 01  LBL "EE"  02  1  03  X<>Y  04  RDN  05  P-R  06  R^  07  X<>Y  08  P-R  09  RDN  10  R-P  11  X<> Z  12  ST- Y  13  X<> Z  14  P-R  15  R^  16  R-P  17  X<>Y  18  RDN  19  R-P  20  X<> T  21  END

( 32 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z e e Y decl b X RA l

where  e = obliquity of the ecliptic

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

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

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

1°)  Gauss-Herrick-Gibbs Method

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

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

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

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

•  R01 = 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 = instant of passage in perihelion  ( in days from 2000/01/01  0h )
R20 = q = perihelion distance    ( in Astronomical Units )
R21 = e = eccentricity
R22 =  i = inclination  ( in degrees )
R23 = omega = argument of perihelion ( in degrees )
R24 = OMEGA = longitude of the ascending node ( in degrees )
R25 = n = mean motion  ( in degrees/day )
R26 = p = parameter ( in AU )

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

R29 = Distance Sun-Comet at the 1st instant  = r1
R28 = (Distance Sun-Comet at the 2nd instant)^3 = r2^3
R30 = Distance Sun-Comet at the 3rd instant      = r3

The coordinates of the vector r1 are stored in R41-R42-R43
The coordinates of the vector r2 are stored in R44-R45-R46         The coordinates of the vector V2 are stored in R34-R35-R36
The coordinates of the vector r3 are stored in R47-R48-R49

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

 01  LBL "GAUSS"   02  DEG   03  STO 50    04  15   05  STO 19   06  49   07  STO 20    08  LBL 01   09  RCL IND 19    10  DSE 19   11  RCL IND 19   12  DSE 19   13  RCL IND 19   14  .5   15  -   16  FC? 01   17  ST- X   18  2807410   19  /   20  23.439279   21  -   22  STO 00    23  RDN   24  HR   25  15   26  *   27  X<>Y   28  HR   29  STO 21   30  COS   31  P-R   32  X<>Y   33  RCL 00   34  X<>Y   35  P-R   36  X<>Y   37  X<> 00   38  RCL 21   39  SIN   40  P-R   41  ST+ 00   42  RDN   43  -   44  RCL 00   45  STO IND 20   46  DSE 20   47  RDN   48  STO IND 20   49  DSE 20   50  X<>Y   51  STO IND 20   52  3   53  ST- 19   54  DSE 20   55  DSE 19   56  GTO 01   57  RCL 44   58  RCL 48   59  *   60  RCL 45   61  RCL 47   62  *   63  -   64  STO 30   65  RCL 46   66  RCL 47   67  *   68  RCL 44   69  RCL 49   70  *   71  -   72  STO 29   73  RCL 45   74  RCL 49   75  *   76  RCL 46   77  RCL 48   78  *   79  -   80  STO 28   81  RCL 42   82  RCL 49   83  *   84  RCL 43   85  RCL 48   86  *   87  -   88  STO 31   89  RCL 43   90  RCL 47   91  *   92  RCL 41   93  RCL 49   94  *   95  -    96  STO 32   97  RCL 41   98  RCL 48   99  * 100  RCL 42 101  RCL 47 102  * 103  - 104  STO 33 105  RCL 42 106  RCL 46 107  * 108  RCL 43  109  RCL 45 110  * 111  - 112  STO 34  113  RCL 47 114  * 115  RCL 43 116  RCL 44  117  * 118  RCL 41 119  RCL 46 120  * 121  - 122  STO 35  123  RCL 48 124  * 125  + 126  RCL 41 127  RCL 45 128  * 129  RCL 42 130  RCL 44 131  * 132  - 133  STO 36 134  RCL 49 135  * 136  + 137  STO 00 138  RCL 13 139  RCL 07 140  - 141  STO 37 142  X^2 143  RCL 07 144  RCL 01 145  - 146  STO 38 147  X^2 148  RCL 13 149  RCL 01 150  - 151  ST/ 37 152  ST/ 38 153  X^2 154  ST- Z 155  X<>Y 156  - 157  RCL 38 158  * 159  STO 40 160  X<>Y 161  CHS 162  RCL 37 163  * 164  3379.380681 165  STO 27 166  6 167  * 168  ST/ 40 169  / 170  STO 39 171  RCL 04 172  * 173  RCL 16 174  RCL 40 175  * 176  + 177  STO M 178  RCL 05 179  RCL 39 180  * 181  RCL 17 182  RCL 40 183  * 184  + 185  STO N 186  RCL 06 187  RCL 39 188  * 189  RCL 18 190  RCL 40 191  * 192  + 193  STO O 194  24 195  STO 25 196  XEQ 02 197  RCL 04 198  RCL 37 199  * 200  RCL 10 201  - 202  RCL 16 203  RCL 38 204  * 205  + 206  STO M 207  RCL 05 208  RCL 37 209  * 210  RCL 11 211  - 212  RCL 17  213  RCL 38  214  * 215  + 216  STO N 217  RCL 06 218  RCL 37 219  * 220  RCL 12  221  - 222  RCL 18 223  RCL 38  224  * 225  + 226  STO O 227  DSE 25 228  XEQ 02 229  " D=?" 230  ASTO O 231  RCL 37 232  RCL 38 233  * 234  RCL 13 235  RCL 01 236  - 237  X^2 238  ST* Y 239  + 240  12 241  RCL 27 242  * 243  STO 26 244  / 245  STO 25 246  RCL 10 247  RCL 44 248  * 249  RCL 11 250  RCL 45 251  * 252  + 253  RCL 12 254  RCL 46 255  * 256  + 257  ST+ X 258  STO M 259  RCL 10 260  X^2 261  RCL 11 262  X^2 263  + 264  RCL 12 265  X^2 266  + 267  STO N 268  GTO 10 269  LBL 02 270  RCL 34 271  RCL M 272  * 273  RCL 35 274  RCL N 275  * 276  + 277  RCL 36 278  RCL O 279  * 280  + 281  RCL 00 282  / 283  STO IND 25 284  DSE 25 285  RCL 31 286  RCL M 287  * 288  RCL 32 289  RCL N 290  * 291  + 292  RCL 33 293  RCL O 294  * 295  + 296  RCL 00 297  / 298  STO IND 25 299  DSE 25 300  RCL 28 301  RCL M 302  * 303  RCL 29 304  RCL N 305  * 306  + 307  RCL 30 308  RCL O 309  * 310  + 311  RCL 00 312  / 313  STO IND 25  314  RTN 315  LBL 10 316  VIEW 50 317  RCL 25  318  RCL 50  319  ENTER 320  X^2 321  * 322  STO 28 323  ST+ Y 324  X^2 325  / 326  STO 29 327  RCL 23  328  * 329  RCL 20 330  + 331  ENTER 332  STO 00 333  RCL M 334  - 335  * 336  RCL N 337  + 338  RCL 50 339  X^2 340  - 341  RCL 25 342  ST+ X 343  RCL 28 344  ST+ Y 345  X^2 346  RCL 50 347  * 348  / 349  RCL 23 350  * 351  3 352  * 353  RCL 00 354  ST+ X 355  RCL M 356  - 357  * 358  RCL 50 359  ST+ X 360  + 361  / 362  ST+ 50 363  RCL 50 364  / 365  ABS 366   E-8 367  XY? 373  X<>Y 374  X>0? 375  GTO 00 376  RCL O 377  TONE 9 378  STOP 379  STO 50 380  GTO 10 381  LBL 00 382  RCL 00 383  ST* 44 384  ST* 45 385  ST* 46 386  RCL 10 387  ST- 44 388  RCL 11 389  ST- 45 390  RCL 12 391  ST- 46 392  RCL 22 393  RCL 29 394  * 395  RCL 19 396  + 397  RCL 39 398  RCL 29 399 * 400   RCL 37 401  + 402  / 403  ST* 41 404  ST* 42 405  ST* 43 406  RCL 04 407  ST- 41 408  RCL 05 409  ST- 42 410  RCL 06 411  ST- 43 412  RCL 24 413  RCL 29 414  * 415  RCL 21 416  + 417  RCL 40 418  RCL 29  419  * 420  RCL 38          421  + 422  / 423  ST* 47 424  ST* 48 425  ST* 49 426  RCL 16  427  ST- 47 428  RCL 17 429  ST- 48 430  RCL 18  431  ST- 49 432  RCL 41 433  X^2 434  RCL 42 435  X^2 436  RCL 43 437  X^2 438  + 439  + 440  SQRT 441  STO 29 442  RCL 47 443  X^2 444  RCL 48 445  X^2 446  RCL 49 447  X^2 448  + 449  + 450  SQRT 451  STO 30 452  RCL 37 453  RCL 07 454  RCL 01 455  - 456  STO 35 457  / 458  STO 31 459  RCL 13 460  RCL 07 461  - 462  ST/ 38 463  RCL 26 464  ST/ 35 465  / 466  STO 34 467  RCL 29 468  3 469  Y^X 470  / 471  + 472  CHS 473  STO 37 474  RCL 41 475  * 476  RCL 35 477  ST- 34 478  RCL 30 479  3 480  Y^X 481  / 482  RCL 38 483  ST- 31 484  + 485  STO 39 486  RCL 47 487  * 488  + 489  RCL 31 490  RCL 34 491  RCL 28 492  / 493  + 494  STO 38 495  RCL 44 496  * 497  + 498  STO 34 499  RCL 42 500  RCL 37 501  * 502  RCL 45 503  RCL 38 504  * 505  + 506  RCL 48 507  RCL 39 508  * 509  + 510  STO 35 511  RCL 43 512  RCL 37 513  * 514  RCL 46 515  RCL 38 516  * 517  + 518  RCL 49 519  RCL 39 520  * 521  + 522  STO 36 523  RCL 45  524  * 525  RCL 46  526  RCL 35          527  * 528  - 529  STO 37  530  X^2 531  RCL 44 532  RCL 36 533  * 534  RCL 46  535  RCL 34 536  * 537  - 538  STO 38 539  X^2 540  + 541  RCL 44 542  RCL 35 543  * 544  RCL 45 545  RCL 34 546  * 547  - 548  STO 39 549  X^2 550  + 551  RCL 27 552  * 553  STO 26 554  RCL 37 555  RCL 38 556  R-P 557  X<>Y 558  STO 24 559  CLX 560  RCL 39 561  R-P 562  X<>Y 563  STO 22 564  RCL 43 565  X<>Y 566  SIN 567  / 568  RCL 41 569  RCL 24 570  COS 571  * 572  RCL 42 573  RCL 24 574  SIN 575  * 576  + 577  R-P 578  X<>Y 579  STO 23 580  RCL 49 581  RCL 22 582  SIN 583  / 584  RCL 47 585  RCL 24 586  COS 587  * 588  RCL 48 589  RCL 24 590  SIN 591  * 592  + 593  R-P 594  CLX 595  RCL 23 596  - 597  2 598  / 599  STO 25 600  RCL 26 601  RCL 29 602  / 603  STO Y 604  RCL 26 605  RCL 30 606  / 607  ST+ Z 608  - 609  RCL 25 610  SIN 611  / 612  X<>Y 613  2 614  - 615  RCL 25 616  COS 617  / 618  R-P 619  2 620  / 621  STO 21 622  CLX 623  RCL 25 624  - 625  ST- 23 626  STO 25 627  RCL 26  628  1 629  RCL 21           630  + 631  / 632  STO 20  633  1 634  RCL 21  635  - 636  X#0? 637  GTO 00 638  X<> 25 639  2 640  / 641  TAN 642  ENTER 643  X^2 644  3 645  + 646  * 647  RCL 20 648  3 649  Y^X 650  ST+ X 651  9 652  / 653  RCL 27 654  * 655  SQRT 656  * 657  GTO 05 658  LBL 00 659  CF 00 660  X<0? 661  SF 00 662  STO Z 663  / 664  3 665  Y^X 666  RCL 27 667  * 668  SIGN 669  LASTX 670  ABS 671  SQRT 672  * 673  1/X 674  R-D 675  X<> 25 676  2 677  / 678  X<>Y 679  1 680  RCL 21 681  + 682  / 683  ABS 684  SQRT 685  FS?C 00 686  GTO 03 687  P-R 688  LASTX 689  / 690  R-P 691  X<>Y 692  ST+ X 693  D-R 694  LASTX 695  SIN 696  GTO 04 697  LBL 03 698  X<>Y 699  TAN 700  * 701  1 702  RCL Y 703  - 704  / 705  ST+ X 706  LN1+X 707  ENTER 708  E^X-1 709  LASTX 710  CHS 711  E^X-1 712  - 713  2 714  / 715  LBL 04 716  RCL 21 717  * 718  - 719  R-D 720  RCL 25 721  / 722  LBL 05 723  RCL 01 724  X<>Y 725  - 726  STO 19 727  CLA 728  END

( 1102 bytes / SIZE 051 )

 STACK INPUT OUTPUT X D T

Where  D  is an estimtion of the distance Sun-Comet at the 2nd instant

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

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

t                 5510                                   5519                                5529                               days
RA        1h07m43s058                     0h58m40s151                   0h53m53s415
Decl       -57°17'23"42                      -52°05'21"91                    -46°54'15"67
X          0.653892160                      0.763553245                   0.863088915                  Astronomical Units
Y        -0.736974521                     -0.624900515                 -0.482202751                   Astronomical Units
Z          0.000019390                      0.000019018                   0.000014378                   Astronomical Units

-Store these 18 numbers into

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

-If you choose  D = 3 AU

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

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

Notes:

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

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

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

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

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

2°)  Position & Velocity >>> Orbital Elements

-If we know the position and velocity vectors r & V  at an instant t , we can deduce the 6 orbital elements  T  q  e  i  omega  OMEGA
-The formulas are given in "Comets for the HP-41"

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

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

>>>> When the program stops,

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

Flag:  F00
Subroutines: /

-Line 100  is a three-byte GTO 02

 01  LBL "RV-EL"   02  DEG   03  STO 00   04  RCL 36    05  RCL 45   06  *   07  RCL 35    08  RCL 46   09  *   10  -   11  STO 37   12  X^2   13  RCL 36   14  RCL 44    15  *   16  RCL 34   17  RCL 46   18  *   19  -   20  STO 38   21  X^2   22  +   23  RCL 35   24  RCL 44   25  *   26  RCL 34   27  RCL 45   28  *   29  -   30  STO 39   31  X^2   32  +   33  3379.380681   34  STO 40    35  *   36  STO 26 37  RCL 37   38  RCL 38   39  R-P   40  X<>Y   41  STO 24    42  CLX   43  RCL 39    44  R-P   45  X<>Y   46  STO 22   47  RCL 46   48  X<>Y   49  SIN   50  /   51  RCL 44          52  RCL 24   53  COS   54  *   55  RCL 45   56  RCL 24   57  SIN   58  *   59  +   60  R-P   61  X<>Y   62  STO 23   63  RCL 34   64  RCL 44   65  *   66  RCL 35   67  RCL 45    68  *   69  +   70  RCL 36   71  RCL 46   72  * 73  +   74  SIGN   75  STO 48   76  2   77  RCL 44    78  X^2   79  RCL 45    80  X^2   81  +   82  RCL 46   83  X^2   84  +   85  SQRT   86  STO 47          87  /   88  RCL 34   89  X^2   90  RCL 35   91  X^2   92  +   93  RCL 36   94  X^2   95  +   96  RCL 40    97  *   98  -   99  X=0? 100  GTO 02 101  CF 00 102  X>0? 103  SF 00 104  1/X 105  STO 27 106  RCL 26 107  LASTX 108  * 109  1 110  - 111  ABS 112  SQRT 113  STO 21  114  1 115  ST+ Y 116  RCL 21 117  - 118  / 119  ABS 120  SQRT 121  STO 25  122  RCL 48        123  RCL 27  124  RCL 47 125  - 126  RCL 27 127  RCL 21 128  * 129  / 130  FS? 00 131  ACOS 132  FS? 00 133  GTO 00 134  ENTER 135  X^2 136  1 137  - 138  SQRT 139  + 140  LN 141  LBL 00 142  * 143  STO 28  144  FS? 00 145  GTO 00 146  E^X-1 147  STO Y 148  2 149  + 150  / 151  RCL 25  152  * 153  ATAN 154  GTO 01 155  LBL 00 156  2 157  / 158  RCL 25        159  P-R 160  LASTX 161  / 162  R-P 163  X<>Y 164  LBL 01 165  ST+ X 166  ST- 23 167  RCL 27 168  3 169  Y^X 170  RCL 40 171  * 172  SIGN 173  LASTX 174  ABS 175  SQRT 176  * 177  1/X 178  R-D 179  STO 25  180  RCL 26 181  1 182  RCL 21 183  + 184  / 185  STO 20 186  RCL 28  187  FS?C 00 188  GTO 00 189  ENTER 190  E^X-1 191  LASTX 192  CHS 193  E^X-1 194  - 195  2 196  / 197  GTO 01 198  LBL 02 199  STO 25        200  RCL 47 201  RCL 26 202  2 203  / 204  STO 20  205  / 206  1 207  STO 21 208  - 209  SQRT 210  RCL 48 211  * 212  ATAN 213  ST+ X 214  ST- 23 215  LASTX 216  ENTER 217  X^2 218  3 219  + 220  * 221  RCL 20  222  3 223  Y^X 224  ST+ X 225  9 226  / 227  RCL 40        228  * 229 SQRT 230  * 231  GTO 03 232  LBL 00 233  D-R 234   LASTX 235  SIN 236  LBL 01 237  RCL 21  238  * 239  - 240  R-D 241  RCL 25  242  / 243  LBL 03 244  RCL 00 245  X<>Y 246  - 247  STO 19 248  END

( 363 bytes / SIZE 049 )

 STACK INPUT OUTPUT X t T

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

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

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

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

and in registers R20 thru R27

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

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

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

and in registers R20 thru R27

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

Notes:

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

3°)  Position & Velocity >>> Position2

-Here, we could use "RV-EL" first to get the orbital elements and then the programs listed in "Comets for the HP-41" to get the position at another instant.
-But, in order to avoid a possible loss of accuracy if the eccentricity is very close to 1, "RV-R" employs the universal variable x  to perform the calculations
( cf reference [2] )

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

here,  z = x2 / a    where   a = semi-major axis              1 / a = 2 / r - V2 / k2

t = ( t2 - t1 )  and  k = 0.01720209895 = Gauss gravitational constant

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

-After founding  x  ,  we have

f = 1 - x2 C(z) / r
g = t - x3 S(z) / k

and   r' = f r + g V

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

•  R44-R45-R46 = r1  = rectangular ecliptic coordinates of the asteroid
•  R34-R35-R36 = V1 = -------------------------------------  velocity    at the same instant t1

R01 thru R06 & R12 thru R26 are unused

>>>> When the program stops,

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

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

-Line 183 is a three-byte  GTO 01

 01  LBL "RV-R"   02  RAD   03  STO 00   04  RCL 34   05  RCL 44           06  *   07  RCL 35   08  RCL 45    09  *   10  +   11  RCL 36   12  RCL 46   13  *   14  +   15  STO 08   16  2   17  RCL 44   18  X^2   19  RCL 45   20  X^2   21  +   22  RCL 46   23  X^2   24  +   25  SQRT   26  STO 07   27   /   28  RCL 34   29  X^2   30  RCL 35   31  X^2   32  +   33  RCL 36 34  X^2   35  +   36  58.13244087   37  STO 09           38  X^2   39  *   40  -   41  STO 27    42  CF 00   43  X<0?   44  SF 00   45  ABS   46  RCL 00   47  *   48  RCL 09   49  ST* 08   50  /   51  STO 10   52  LBL 01   53  RCL 10   54  X^2   55  RCL 27   56  *   57  STO 11   58  ABS   59  .73   60  X<=Y?   61  GTO 03   62  4   63  ENTER   64  CLX   65  LBL 02   66  RCL Y 67  ST+ X   68  3   69  +   70  FACT   71  1/X   72  +   73  RCL 11           74  *   75  CHS   76  DSE Y   77  GTO 02   78  6   79  1/X   80  +   81  GTO 04   82  LBL 03   83  X<>Y   84  SQRT   85  ENTER   86  ENTER   87  FS? 00   88  GTO 03   89  SIN   90  GTO 00   91  LBL 03   92  E^X-1   93  LASTX   94  CHS   95  E^X-1   96  -   97  2   98  /   99  LBL 00 100  - 101  X<>Y 102  / 103  RCL 11         104  / 105  LBL 04 106  RCL 11  107  ABS 108  SQRT 109  2 110  / 111  FS? 00 112  GTO 03 113  SIN 114  X^2 115  ST+ X 116  RCL 11 117  GTO 04 118  LBL 03 119  E^X-1 120  LASTX 121  CHS 122  E^X-1 123  - 124  X^2 125  RCL 11 126  ST+ X 127  CHS 128  LBL 04 129  X#0? 130  GTO 04 131  RDN 132  SIGN 133  2 134  LBL 04 135  / 136  STO 28         137  RCL 08 138  * 139  X<>Y 140  STO 29 141  RCL 10  142  * 143  + 144  RCL 10 145  * 146  1 147  RCL 11 148  RCL 29 149  * 150  - 151  STO Z 152  RCL 07 153  * 154  + 155  RCL 10 156  * 157  RCL 00 158  RCL 09 159  / 160  - 161  RCL 10 162  RCL 28 163  * 164  R^ 165  RCL 08 166  * 167  + 168  RCL 10         169  * 170  RCL 11  171  RCL 28 172  * 173  1 174  - 175  RCL 07  176  * 177  - 178  / 179  ST- 10 180  ABS 181   E-7 182  X 29 203  RCL 10         204  3 205  Y^X 206  * 207  RCL 09  208  * 209  RCL 00 210  X<>Y 211  - 212  RCL 34 213  X<>Y 214  * 215  ST+ 27 216  RCL 35 217  LASTX 218  * 219  ST+ 28 220  RCL 36 221  LASTX 222  * 223  ST+ 29 224  RCL 29 225  RCL 28 226  RCL 27 227  DEG 228  CF 00 229  END

( 310 bytes / SIZE 047 )

 STACK INPUT OUTPUT Z / Z Y / Y X t2 - t1 X

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

Example1:        r1  = ( 0.16 , 1.38 , 0.24 )    V1 = ( 0.015 , 0.01 , 0.001 )    Find the position 100 days later

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

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

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

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

Notes:

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

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

 44  GTO 00   45  RCL 00   46  *   47  RCL 09   48  ST* 08   49  /   50  STO 10   51  GTO 01   52  LBL 00 53  SF 00   54  RCL 00   55  ST+ X   56  *   57  CHS   58  1   59  RCL 07   60  RCL 27   61  * 62  -   63  RCL 27   64  CHS   65  SQRT   66  RCL 09   67  *   68  /   69  RCL 00   70  SIGN 71  *   72  RCL 08   73  +   74  /   75  RCL 09    76  ST* 08   77  X^2   78  /   79  X#0? 80  LN   81  RCL 27   82  CHS   83  SQRT   84  /   85  RCL 00    86  SIGN   87  *   88  STO 10

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

4°)  2 Positions >>> Velocity

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

-You have to give an estimate of the parameter p

-Flag F04 must be cleared if the angle  ( r1 ; r2 ) < 180°
-Flag F04 must be set if the angle  ( r1 ; r2 ) > 180°

-Let    C = Cos ( v2 - v1 ) = ( r1.r2 ) / ( r1. r2 )          where    v  is the true anomaly                  Sin ( v2 - v1 ) = +/- ( 1 - C2 )1/2

-We have   f = ( r2 / p ) ( C - 1 ) + 1   &  g = ( r1. r2 ) / ( k p1/2 )  Sin ( v2 - v1 )

from which   V1 =  ( r2 - f r1 ) / g   whence  an estimate of r2 ( say r* )   whence   C* = ( r1.r* ) / ( r1. r* )

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

-This equation is solved by the secant method.

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

•  R44-R45-R46 = r1  = rectangular ecliptic coordinates of the asteroid at the instant t1
•  R41-R42-R43 = r2 = ----------------------------------------------- at the instant t2

>>>> When the program stops,   R34-R35-R36 = V1 = rectangular ecliptic coordinates of the velocity at the first instant t1

Flag:  F04    Clear F04 if the angle  ( r1 ; r2 ) < 180°
Set F04 if the angle  ( r1 ; r2 ) > 180°

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

 01  LBL "RR-V"   02  STO 03    03  .99   04  /   05  STO 04    06  RDN   07  STO 02    08  X<>Y   09  STO 01         10  RCL 41    11  RCL 44           12  *   13  RCL 42   14  RCL 45   15  *   16  +   17  RCL 43   18  RCL 46   19  *   20  +   21  RCL 44   22  X^2   23  RCL 45   24  X^2   25  +   26  RCL 46   27  X^2   28  +   29  SQRT   30  STO 31 31  RCL 41   32  X^2   33  RCL 42    34  X^2   35  +   36  RCL 43    37  X^2   38  +   39  SQRT   40  STO 32         41  *   42  /   43  STO 30   44  XEQ 10   45  STO 33           46  RCL 03   47  X<> 04   48  STO 03   49  LBL 01   50  XEQ 10   51  ENTER   52  ENTER   53  X<> 33   54  -   55  X#0?   56  /   57  RCL 04   58  RCL 03   59  STO 04   60  STO T 61  -   62  *   63  +   64  STO 03    65  LASTX   66  X<>Y   67  /   68  ABS   69  VIEW 03   70   E-7   71  XY   86  *   87  STO 34   88  RCL 45   89  LASTX   90  * 91  STO 35   92  RCL 46   93  LASTX   94  *   95  STO 36    96  RCL 41   97  ST+ 34   98  RCL 42    99  ST+ 35 100  RCL 43 101  ST+ 36 102  RCL 30         103  ACOS 104  SIN 105  FS? 04 106  CHS 107  RCL 31 108  RCL 32 109  * 110  * 111  RCL 03  112  SQRT 113  / 114  58.13244087 115  * 116  ST/ 34 117  ST/ 35 118  ST/ 36 119  RCL 02 120  RCL 01 121  - 122  XEQ "RV-R" 123  RCL 44 124  RCL 27  125  * 126  RCL 45 127  RCL 28  128  * 129  + 130  RCL 46          131  RCL 29        132  * 133  + 134  RCL 31 135  / 136  RCL 27 137  X^2 138  RCL 28 139  X^2 140  + 141  RCL 29 142  X^2 143  + 144  SQRT 145  / 146  RCL 30 147  - 148  RTN 149  END

( 286 bytes / SIZE 049 )

 STACK INPUTS OUTPUTS1 Z t1 / Y t2 / X p0 p

When the program stops,  R03 = p

Example:    t1 = 0   ,   t1 = 100 days  ,    r1 = ( 0.16 , 1.38 , 0.24 )  ,  r2 = ( 1.509299637 , 1.919542031 , 0.265117223 )

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

-If you choose  p0 = 1   and   CF 04

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

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

Notes:

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

-Try another p-estimation ( often a smaller one ) and execute "RR-V" again.

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

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

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

-Assuming there are n observations, you have to give an estimate of the distances Earth-Asteroid at the first and last instants  D1 & Dn
and an estimate of the parameter p.
-So, the HP-41 can compute the position-vectors ri  for  i = 2 , 3 , .... , n-1  thanks to "RR-V" & "RV-R"

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

Pi = ( ri + Ri ).Ai               where      Ri = rectangular ecliptic coordinates of the Sun                                                  Li = eclptic longitude
Qi = ( ri + Ri ).Di                       Ai = ( -Sin Li , Cos Li , 0 ) Di = ( -Sin bi Cos Li , -Sin bi Sin Li , Cos bi )      bi = ecliptic latitude

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

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

(Pi/D1)DD1 + (Pi/Dn)DDn = - Pi
(Qi/D1)DD1 + (Qi/Dn)DDn = - Qi

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

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

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

•  R51 = t1               •  R57 = t2         ...............    •  R45+6n = tn                  expressed in days since 2000/01/01  0h TT
•  R52 = L1              •  R58 = L2        ...............    •  R46+6n = Ln                Ecliptic longitudes in degrees
•  R53 = b1              •  R59 = b2         ...............   •  R47+6n = bn                 Ecliptic latitudes in degrees
•  R54 = X1             •  R60 = X2         ..............    •  R48+6n = Xn                 in Astronomical Units          rectangular ecliptic
•  R55 = Y1             •  R61 = Y2         ..............    •  R49+6n = Yn                 ---------------------              coordinates
•  R56 = Z1             •  R62 = Z2          ..............    •  R50+6n = Zn                 ---------------------               of the Sun

>>>> When the program eventually stops,

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

Flags:  F00 - F04      Clear F04 if the angle  ( r1 ; rn ) < 180°
Set F04 if the angle  ( r1 ; rn ) > 180°

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

-Line 148 is a three-byte  GTO 10

 01  LBL "HRG"   02  STO 06   03  RDN   04  STO 05    05  RDN   06  STO 03   07  X<>Y   08  STO 21   09  LBL 10   10  RCL 50             11  FRC   12   E3   13  *   14  INT   15  1   16  +   17  STO 24    18  XEQ 01   19  RCL 21    20  ST+ 05   21  XEQ 01   22  RCL 21    23  ST- 05   24  ST+ 06   25  XEQ 01   26  RCL 21   27  ST- 06   28  CLX   29  STO 14   30  STO 15   31  STO 16   32  STO 17   33  STO 18   34  ISG 50   35  .006   36  ST- 50   37  87.041006   38  STO 24   39  LBL 02   40  RCL 24   41  REGMOVE   42  XEQ 03   43  STO 12   44  X<>Y   45  STO 22   46  RCL 24   47  6 48  +   49  REGMOVE   50  XEQ 03   51  RCL 12    52  -   53  RCL 21   54  /   55  STO 13   56  RCL 12   57  *   58  ST+ 14   59  X<>Y   60  RCL 22             61  -   62  RCL 21   63  /   64  STO 23    65  RCL 22    66  *   67  ST+ 14   68  RCL 13   69  X^2   70  RCL 23    71  X^2   72  +   73  ST+ 16   74  RCL 24   75  12   76  +   77  REGMOVE   78  XEQ 03   79  RCL 12   80  -   81  RCL 21   82  /   83  X^2   84  ST+ 18   85  LASTX   86  ST* 13   87  RCL 12   88  *   89  ST+ 15   90  R^   91  RCL 22   92  -   93  RCL 21   94  / 95  X^2   96  ST+ 18   97  LASTX   98  ST* 23   99  RCL 22  100  * 101  ST+ 15 102  RCL 23 103  RCL 13 104  + 105  ST+ 17 106  ISG 50 107  GTO 02 108  RCL 50           109  .006 110  + 111  FRC 112  51 113  + 114  STO 50  115  RCL 15 116  RCL 17  117  * 118  RCL 14  119  RCL 18 120  * 121  - 122  RCL 16 123  RCL 18 124  * 125  RCL 17 126  X^2 127  - 128  STO 12 129  / 130  STO 13 131  ST+ 05 132  RCL 14 133  RCL 17 134  * 135  RCL 15 136  RCL 16 137  * 138  - 139  RCL 12 140  / 141  STO 14 142  ST+ 06 143  RCL 03 144  RCL 06 145  RCL 05  146  RTN 147  X#0? 148  GTO 10 149  XEQ 01 150  RCL 51 151  XEQ "RV-EL" 152  RTN 153  LBL 01 154  RCL 53           155  RCL 52 156  RCL 05 157  XEQ "S-R" 158  RCL 54 159  - 160  STO 44  161  CLX 162  RCL 55  163  - 164  STO 45 165  CLX 166  RCL 56  167  - 168  STO 46 169  RCL 50 170  FRC 171   E3 172  * 173  INT 174  3 175  - 176  STO 04 177  RCL IND X 178  DSE Y 179  RCL IND Y 180  RCL 06 181  XEQ "S-R" 182  ISG 04 183  CLX 184  RCL IND 04 185  - 186  STO 41 187  ISG 04 188  CLX 189  CLX 190  RCL IND 04 191  - 192  STO 42  193  ISG 04 194  CLX 195  CLX 196  RCL IND 04 197  - 198  STO 43 199  5 200  ST- 04 201  RCL 51           202  RCL IND 04 203  RCL 03 204  XEQ "RR-V" 205  CLD 206  RCL 34  207  STO 41 208  RCL 35 209  STO 42  210  RCL 36 211  STO 43 212  RCL 24  213  INT 214  .006 215  + 216   E3 217  / 218  41 219  + 220  REGMOVE 221  6 222  ST+ 24 223  RTN 224  LBL 03 225  RCL 41 226  STO 34 227  RCL 42 228  STO 35 229  RCL 43 230  STO 36 231  RCL IND 50 232  RCL 51 233  - 234  XEQ "RV-R" 235  RCL 50 236  INT 237  5 238  + 239  RCL IND X 240  ST+ 29 241  DSE Y 242  RCL IND Y 243  ST+ 28 244  DSE Z 245  RCL IND Z 246  ST+ 27 247  R^ 248  DSE X 249  RCL IND X 250  1 251  P-R 252  STO 04           253  RDN 254  STO 07 255  DSE Y 256  RCL IND Y 257  LASTX 258  P-R 259  STO 08  260  RCL 27  261  * 262  X<>Y 263  STO 09 264  RCL 28 265  * 266  + 267  RCL 07 268  * 269  RCL 29 270  RCL 04 271  * 272  X<>Y 273  - 274  RCL 28 275  RCL 08 276  * 277  RCL 27 278  RCL 09 279  * 280  - 281  RTN 282  END

( 471 bytes / SIZE 69+6n )

 STACK INPUTS OUTPUTS1 OUTPUTS2 OUTPUTS3 ... OUTPUT T h / / / / Z p / / / / Y D1 D'n D"n D"'n  ... / X Dn D'1 D"1 D"'1   ... T

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

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

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

1st Step          Store these 36 numbers into

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

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

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

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

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

R32 = D1 = 2.316488230

-Likewise, for the last 3 observations

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

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

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

R33 = Dn = 2.717138781

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

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

 01  LBL "INIT"   02  6  03  *  04  50.06  05  +  06   E3  07  /  08  51  09  +  10  STO 49 11  STO 50   12  LBL 01  13  RCL IND 49  14  .5  15  -  16  1  17  ST+ 49  18  CLX  19  RCL IND 49  20  HR 21  15  22  *  23  1  24  ST+49  25  CLX  26  RCL IND 49  27  HR  28  RCL Z  29  2807410  30  / 31  23.439279  32  X<>Y  33  FC? 01  34  CLX  35  -  36  X<> Z  37  XEQ "EE"  38  X<>Y  39  STO IND 49  40  1 41  ST- 49  42  X<> Z  43  STO IND 49  44  1  45  ST- 49  46  ISG 49  47  GTO 01  48  RCL 50  49  END

( 101 bytes )

 STACK INPUT OUTPUT X n 51.eee06

where n = number of observations

Example:     Here  n = 6  so,

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

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

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

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

CF 04

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

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

R/S    >>>>   D1 = 2.314978136 = R05
X<>Y   Dn = 2.715012525 = R06

R/S    >>>>   D1 = 2.314977837 = R05
X<>Y   Dn = 2.715012172 = R06

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

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

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

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

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

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

-The JPL gives almost the same results:

Element                  Value                                    Uncertainty (1-sigma)    Units

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

Notes:

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

6°)  A Test

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

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

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

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

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

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

•  R51 = t1               •  R57 = t2         ...............    •  R45+6n = tn                  expressed in days since 2000/01/01  0h TT
•  R52 = L1              •  R58 = L2        ...............    •  R46+6n = Ln                Ecliptic longitudes in degrees
•  R53 = b1              •  R59 = b2         ...............   •  R47+6n = bn                 Ecliptic latitudes in degrees
•  R54 = X1             •  R60 = X2         ..............    •  R48+6n = Xn                 in Astronomical Units          rectangular ecliptic
•  R55 = Y1             •  R61 = Y2         ..............    •  R49+6n = Yn                 ---------------------              coordinates
•  R56 = Z1             •  R62 = Z2          ..............    •  R50+6n = Zn                 ---------------------               of the Sun

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

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

 01  LBL "EPS"  02  CLX  03  STO 49  04  RCL 50  05  STO 48  06  LBL 01  07  RCL 51  08  FS? 02  09  RCL 57   10  RCL IND 48  11  X<>Y 12  -  13  XEQ "RV-R"  14  RCL 48   15  INT  16  5  17  +  18  RCL IND X   19  ST+ 29  20  DSE Y  21  RCL IND Y  22  ST+ 28 23  DSE Z  24  RCL IND Z   25  ST+ 27  26  RCL 29   27  RCL 28  28  RCL 27  29  XEQ "R-S"  30  CLX  31  2  32  RCL 48  33  INT 34  +  35  X<> Z   36  RCL IND Z   37  ST- Y  38  COS  39  DSE T  40  RCL IND T   41  R^  42  -  43  SIN  44  ASIN 45  *  46  X^2  47  X<>Y  48  X^2  49  +  50  ST+ 49  51  ISG 48  52  GTO 01  53  RCL 49         54  END

( 100 bytes / SIZE 51+6n )

 STACK INPUT OUTPUT X / eps

where  eps = SUM [ dL2Cos2b + db2 ]

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

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

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

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

Notes:

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

RCL 50   FRC  ISG X   INT  /   SQRT

-The formula  dL2Cos2b + db2 for the square of the angular separations is only valid if the separation is small.

References: