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

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

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	/	(x²+y²)^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     l = 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	/	(x²+y²)^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 b cos l = cos d cos a
             cos b sin l = sin e sin d + cos e cos d sin a

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

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 V₂ is computed too.

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

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

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

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

>>>> When the program stops,

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

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

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

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

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

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 X<Y?
368 GTO 10
369 CLD
370 RCL 00
371 RCL 50
372 X>Y?
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) = x³ S(z) + ( r.V / k ) x² C(z) + r x ( 1 - z S(z) ) - t k = 0 by the method of Newton

here, z = x² / a where a = semi-major axis 1 / a = 2 / r - V² / k²

t = ( t₂ - t₁ ) and k = 0.01720209895 = Gauss gravitational constant

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

-After founding x , we have

f = 1 - x² C(z) / r
g = t - x³ S(z) / k

and r' = f r + g V

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

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

R01 thru R06 & R12 thru R26 are unused

>>>> When the program stops,

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

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

-Line 183 is a three-byte GTO 01

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<Y?
183 GTO 01
184 RCL 44
185 1
186 RCL 10
187 X^2
188 RCL 28
189 *
190 RCL 07
191 /
192 -
193 *
194 STO 27
195 RCL 45
196 LASTX
197 *
198 STO 28

199 RCL 46
200 LASTX
201 *
202 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	t₂ - t₁	X

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

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

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

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

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

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

Notes:

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

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

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 r₁ & r₂ at different instants t₁ & t₂
"RR-V" returns the velocity vector V₁ at the first instant
-It uses an iteration method called the method of Herrick & Liu or p-iteration method ( cf reference [3] )

-You have to give an estimate of the parameter p

-Flag F04 must be cleared if the angle ( r₁ ; r₂ ) < 180°
-Flag F04 must be set if the angle ( r₁ ; r₂ ) > 180°

-Let C = Cos ( v₂ - v₁ ) = ( r₁.r₂ ) / ( r₁. r₂ ) where v is the true anomaly Sin ( v₂ - v₁ ) = +/- ( 1 - C² )^1/2

-We have f = ( r₂ / p ) ( C - 1 ) + 1 & g = ( r₁. r₂ ) / ( k p^1/2 ) Sin ( v₂ - v₁ )

from which V₁ = ( r₂ - f r₁ ) / g whence an estimate of r₂ ( say r* ) whence C* = ( r₁.r* ) / ( r₁. r* )

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

-This equation is solved by the secant method.

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

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

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

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

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

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 X<Y?
72 GTO 01
73 RCL 03
74 RTN
75 LBL 10
76 1
77 RCL 30
78 -
79 RCL 32
80 *
81 RCL 03
82 ST- Y
83 /
84 RCL 44
85 X<>Y
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	t₁	/
Y	t₂	/
X	p₀	p

When the program stops, R03 = p

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

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

-If you choose p₀ = 1 and CF 04

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

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

Notes:

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

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

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

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

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

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

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

P_i = ( r_i + R_i ).A_i where R_i = rectangular ecliptic coordinates of the Sun L_i = eclptic longitude
Q_i = ( r_i + R_i ).D_i A_i = ( -Sin L_i , Cos L_i , 0 ) D_i = ( -Sin b_i Cos L_i , -Sin b_i Sin L_i , Cos b_i ) b_i = ecliptic latitude

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

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

(¶P_i/¶D₁)DD₁ + (¶P_i/¶D_n)DD_n = - P_i
(¶Q_i/¶D₁)DD₁ + (¶Q_i/¶D_n)DD_n = - Q_i

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

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

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

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

>>>> When the program eventually stops,

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

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

-Line 148 is a three-byte GTO 10

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	D₁	D'_n	D"_n	D"'_n ...	/
X	D_n	D'₁	D"₁	D"'₁ ...	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 = D₁ = 2.316488230

-Likewise, for the last 3 observations

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

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

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

R33 = D_n = 2.717138781

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

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

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 D₁ = 2.316488230 D_n = 2.717138781 thus,

CF 04

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

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

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

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

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

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

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

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

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

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

-The JPL gives almost the same results:

Element Value Uncertainty (1-sigma) Units

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

Notes:

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

6°) A Test

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

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

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

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

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

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

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 [ dL²Cos²b + db² ]

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

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

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

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

Notes:

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

RCL 50 FRC ISG X INT / SQRT

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

References:

[1] "SOLEX" which may be downloaded from http://chemistry.unina.it/~alvitagl/solex/
[2] "FUNDAMENTALS OF ASTRODYNAMICS" - Roger R. Bate Donald D. Mueller Jerry E. White
[3] "Fundamentals of Celestial Mechanics" - JMA Danby - Willmann-Bell - ISBN 0-943396-20-4
[3] Robin M. Green - "Spherical Astronomy" - Cambridge University Press - ISBN 0-521-31779-7
[4] John D. Anderson - JPL Technical Report n° 32-497
[5] Dan Boulet - "Methods of Orbit Determination for the Micro Computer" - Willmann-Bell - ISBN 0-943396-34-4

hp41programs

Orbit Determination for the HP-41