hp41programs

Kepler

Kepler's Equation for the HP-41


Overview
 

 1°)  A Simple Method

   a)   Focal Program
   b)   M-Code Routine

 2°)  Newton's Method

   a)   Focal Program
   b)   M-Code Routine

 3°)  Better rate of convergence

   a)   Focal Program
   b)   M-Code Routine

 4°)  Approximate Formula ( e < 1 )
 

-The following programs solve the equation of Kepler:

      M = E - e sin E     if  e < 1       ( elliptic case )
      M = e sinh E - E   if  e > 1       ( hyperbolic case )

 where  e = eccentricity ,  M = mean anomaly & E = eccentric anomaly

-If e = 1 ( parabolic case ) cf "Comets for the HP-41"
 

1°)  A  Simple Method
 

     a) Focal Program
 

-The simplest method to solve Kepler's equation is the iteration                                            En+1 = M + e sin En                    starting with    E0 = M      if  e < 1
-In the hyperbolic case, the equation is re-written   E = Arc sinh [ ( M+E ) / e ]  whence      En+1 = Arc sinh [ ( M+En ) / e ]     also with       E0 = M      if  e > 1

-Though the convergence is slow, it works well provided  e is not too close to 1, for example the 9 major planets.
 

Data Registers:    R00 = M ,  R01 = e or  180 e / PI
Flags: /
Subroutines: /
 
 

01  LBL "KPL"
02  STO 00       
03  CLX
04  SIGN
05  X<>Y
06  X>Y?
07  GTO 00
08  DEG
09  R-D
10  STO 01
11  RCL 00       
12  LBL 01
13  ENTER^
14  SIN
15  RCL 01
16  *
17  RCL 00
18  +
19  X#Y?
20  GTO 01 
21  RTN
22  LBL 00       
23  STO 01
24  RCL 00
25  LBL 02
26  STO Y
27  RCL 00
28  +
29  RCL 01       
30  /
31  ENTER^
32  X^2
33  SIGN
34  LASTX
35  +
36  SQRT
37  +
38  LN
39  X#Y?
40  GTO 02       
41  END

 
  ( 53 bytes / SIZE 002 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             e            E
           X            M            E

 
Examples:

  e = 0.25   ENTER^
 M = 12°    XEQ "KPL"  >>>>  E = 15°93182894         ---Execution time = 27s---

  e = 1.2     ENTER^
 M =  3      XEQ "KPL"   >>>>  E = 2.166183261         ---Execution time = 11s---

Notes:

- M is expressed in degrees only if e < 1.

-Infinite loops are possible because in this program, the iterations stop when 2 consecutive approximations are equal.
-For instance, if e = 0.999999  and  M = 179°99  the successive results are  179.99  180  179.99  180  and so on.
-Moreover, if e is close to 1 and M is small, execution time becomes prohibitive.
-Newton's method is of course better in those cases.
 

     b) M-Code Routine
 

-This routine does not check for alpha data
 

08C  "L"
010   "P"
00B  "K"
2A0   SETDEC
0F8   READ 3(X)
128   WRIT 4(L)
2EE   ?C#0 ALL
3A0   ?NC RTN
068   WRIT 1(Z)
0B8   READ 2(Y)
00E   A=0 ALL            A
1BE  A=A-1 MS
35C  PT=12                 =
162   A=A+1 @PT      -1
01D  ?NCXQ              C=
060   1807                  A+C
104   CLRF 8
2FE   ?C#0 MS
113   JNC +22h=+34d
108   SETF 8
0B8   READ 2(Y)
231   ?NCXQ             performs the
064   198C              R-D  conversion
028   WRIT 0(T)
0F8   READ 3(X)
070   N=C ALL            C
3C4  ST=0
229   ?NCXQ               =
048   128A                sin(C°)
10E   A=C ALL
046   C=0 S&X            C
270   RAMSLCT          =
038   READDATA       T
135   ?NCXQ              C=
060   184D                 A*C
10E   A=C ALL
138   READ 4(L)
01D  ?NCXQ              C=
060   1807                  A+C
0EE   C<>B ALL
078   READ 1(Z)
10E   A=C ALL
0F8   READ 3(X)
068   WRIT 1(Z)
0EE   C<>B ALL
0E8   WRIT 3(X)
36E   ?A#C ALL
10B   JNC +21h=+33d                This line may be simply replaced by   3A0   ?NC RTN
3CC  ?KEY                                  in this case, delete the last 11 instructions written in red
360   ?C RTN
10C   ?FSET 8
337    JC -1Ah=-26d
0F8   READ 3(X)
10E   A=C ALL
138   READ 4(L)
01D  ?NCXQ              C=
060   1807                  A+C
10E   A=C ALL
0B8   READ 2(Y)
261   ?NCXQ              C=
060   1898                  A/C
070   N=C ALL
10E   A=C ALL
135   ?NCXQ              C=
060   184D                 A*C
00E   A=0 ALL            A
35C  PT=12                 =
162   A=A+1 @PT      1
01D  ?NCXQ              C=
060   1807                  A+C
2F9   ?NCXQ              C=
060   18BE                sqrt(C)
10E   A=C ALL
0B0   C=N ALL
01D  ?NCXQ              C=
060   1807                  A+C
084   CLRF 5               C
115   ?NCXQ               =
06C  1B45                  Ln(C)
2C3  JNC -28h=-40d
078   READ 1(Z)
01D  ?NCXQ
060   1807
10E   A=C ALL
04E   C=0 ALL
35C   PT=12
090   LD@PT- 2
261   ?NCXQ
060   1898
0E8   WRIT 3(X)
3E0   RTN
 
 

      STACK        INPUTS      OUTPUTS
           T             /    /   or  180e/PI
           Z             /         temp
           Y             e            e
           X            M            E
           L             /            M

 
Examples:

  e = 1.2     ENTER^
 M =  3      XEQ "KPL"   returns   E = 2.166183261   in 7 seconds

  e  = 1.2     ENTER^
 M = 0.4     XEQ "KPL"  returns   E = 0.987853766   in 18 seconds

  e = 0.25   ENTER^
 M = 12°    XEQ "KPL"  returns   E = 15°93182894    in 9 seconds

  e  = 0.9     ENTER^
 M =  7°      XEQ "KPL"  returns   E = 40°46693450   in 35 seconds

Notes:

-To avoid infinite loops, the routine stops when  En = En-2  and the last 11 lines compute  ( En + En-1 )/2
-Thus,  "KPL"  correctly gives  E = 179°995  if  e = 0.999999  and  M = 179°99

-However, despite these tricks, the program remains very slow when  e is close to 1 and M is small.
-So, Newton's method is by far preferable if you have to deal with near-parabolic orbits.
 

2°)  Newton's Method
 

     a) Focal Program
 

-Kepler's equation is re-written under the form  f(E) = 0  and the successive approximations are    En+1 = En - f( En) / f '(En)  again with E0 = M
 

Data Registers:   R00 = M , R01 = E , R02 = e ,  R03: temp ( only if e > 1 )
Flags: /
Subroutines: /
 
 

01  LBL "KPL"
02  DEG
03  STO 00 
04  STO 01         
05  X<>Y
06  STO 02
07  LBL 01
08  CLX
09  SIGN
10  RCL 02
11  X>Y?
12  GTO 02
13  RCL 01         
14  X<>Y
15  P-R
16  ST- Z
17  RDN
18  R-D
19  RCL 01
20  -
21  RCL 00
22  +
23  X<>Y
24  GTO 00
25  LBL 02
26  RCL 01         
27  RCL 00 
28  RCL 01
29  +
30  R^
31  /
32  ENTER^
33  X^2
34  SIGN
35  LASTX
36  +
37  SQRT
38  STO 03         
39  +
40  LN
41  -
42  RCL 03
43  1/X
44  RCL 02
45  ST* Z
46  -
47  LBL 00         
48  /
49  ST+ 01
50  ABS
51   E-6
52  X<Y?
53  GTO 01
54  RCL 01
55  END

 
  ( 71 bytes / SIZE 003 or 004 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             e            /
           X            M            E

 
Examples:

  e = 0.25   ENTER^
 M = 12°    XEQ "KPL"  >>>>  E = 15°93182894         ---Execution time = 4s---

  e = 1.2     ENTER^
 M =  3      XEQ "KPL"   >>>>  E = 2.166183261         ---Execution time = 5.6s---
 

Notes:

-The iterations stop when the difference between 2 consecutive approximations is smaller than 10 -6  ( line 51 )
-Change this line if a lesser accuracy is acceptable.
-A better starting value may be used: cf "Astronomical Algorithms" or the page "Comets for the HP-41"
 

     b) M-Code Routine
 

-This routine does not check for alpha data
 

08C  "L"
010   "P"
00B  "K"
2A0   SETDEC
0B8   READ 2(Y)
10E   A=C ALL
0F8   READ 3(X)
128   WRIT 4(L)
2EE   ?C#0 ALL
3A0   ?NC RTN
04E    C=0 ALL                 C
2BE   C=-C-1 MS
35C   PT= 12
050    LD@PT- 1               =
21C   PT=2
250    LD@PT- 9
250    LD@PT- 9
110    LD@PT- 4          - 10^(-6)
028    WRIT 0(T)              T=C
046    C=0 S&X          C becomes -1
01D   ?NCXQ                   C=
060    1807                      A+C
104    CLRF 8
2FE   ?C#0 MS
1FB   JNC +3Fh=+63d
108    SETF 8
0B8    READ 2(Y)
10E    A=C ALL
0F8    READ 3(X)
0A8   WRIT 2(Y)
0AE   A<>C ALL
0E8   WRIT 3(X)
070    N=C ALL           performs
3C4   ST=0                      the
1D5   ?NCXQ         polar-rectangular          ( in DEG mode without changing the angular mode )   The results are in CPU registers C & N
078    1E75                  conversion
2BE   C=-C-1 MS        C=-C
00E    A=0 ALL               A
35C   PT=12                    =
162    A=A+1 @PT         1
01D   ?NCXQ                 C=
060    1807                    A+C
068    WRIT 1(Z)
0B0   C=N ALL
231    ?NCXQ             performs the
064    198C              R-D  conversion
0EE    C<>B ALL
0B8    READ 2(Y)
10E    A=C ALL
0F8    READ 3(X)
0A8   WRIT 2(Y)
0AE   A<>C ALL
0E8    WRIT 3(X)
2BE   C=-C-1 MS        C=-C
06E    A<>B ALL
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
138    READ 4(L)
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
078    READ 1(Z)
261    ?NCXQ                C=
060    1898                    A/C
070    N=C ALL
10E    A=C ALL
0F8    READ 3(X)
01D   ?NCXQ                 C=
060    1807                    A+C
0E8    WRIT 3(X)
0B0    C=N ALL
05E    C=0 MS             C = | C |
10E    A=C ALL
046    C=0 S&X               C
270    RAMSLCT             =
038    READDATA          T
01D   ?NCXQ                 C=
060    1807                    A+C
2FE    ?C#0 MS
360    ?C RTN
3CC   ?KEY                   these 2 lines allow you to stop an infinite loop
360    ?C RTN                by pressing any key
10C   ?FSET 8
237    JC -3Ah=-58d
013    JNC +02
33B    JNC -19h=-25d
0F8    READ 3(X)
10E    A=C ALL
138    READ 4(L)
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
0B8    READ 2(Y)
261    ?NCXQ               C=
060    1898                    A/C
070    N= C ALL
10E    A=C ALL
135    ?NCXQ                C=
060    184D                   A*C
00E    A=0 ALL               A
35C   PT=12                    =
162    A=A+1 @PT         1
01D   ?NCXQ                 C=
060    1807                    A+C
2F9    ?NCXQ                C=
060    18BE                 sqrt(C)
268    WRIT 9(Q)
22D   ?NCXQ                C=
060    188B                    1/C
10E    A=C ALL
0B8    READ 2(Y)
2BE   C=-C-1 MS        C=-C
01D   ?CXQ                   C=
061    1807                    A+C
068    WRIT 1(Z)
278    READ 9(Q)
10E    A=C ALL
0B0    C=N ALL
01D   ?NCXQ                 C=
060    1807                    A+C
084    CLRF 5                 C
115    ?NCXQ                 =
06C   1B45                   Ln(C)
2BE   C=-C-1 MS        C=-C
10E    A=C ALL
0F8    READ 3(X)
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
0B8    READ 2(Y)
135    ?NCXQ                C=
060    184D                   A*C
28B    JNC -2Fh=-47d
 
 
 

      STACK        INPUTS      OUTPUTS
           T             /      -1 E-06
           Z             /         temp
           Y             e            e
           X            M            E
           L             /            M

 
Examples:

  e  = 1.2     ENTER^
 M = 0.4     XEQ "KPL"  returns   E = 0.987853767   in 3.5 seconds

  e  = 0.9     ENTER^
 M =  7°      XEQ "KPL"  returns   E = 40.46693453°   in 5 seconds
 

3°)  Better Rate of Convergence
 

     a) Focal Program
 

-Newton's method still presents some difficulties. For example, with e = 0.999 and M = 7° , the M-Code routine above gives the correct result ... in 69 seconds!

-The program hereunder uses the following formula to solve f(E) = 0:

      En+1 = En -  2 f( En) / { f '(En) - [ f '(En)2  - 2 f( En) f ''(En) ] 1/2 }              ( again with E0 = M )

      f(E) = 0  is written so that  f ' < 0

-Newton's method uses a linear approximation of f
-Here, we employ a Taylor series up to the terms of degree 2
 

Data Registers:  R00 = M    R02 = e          R05 = - 2 f( E) f ''(E)
                             R01 = E     R03 = f '(E)
Flags: /
Subroutines: /
 
 

01  LBL "KPL"
02  DEG
03  STO 00
04  STO 01
05  X<>Y
06  STO 02
07  LBL 01         
08  RCL 02
09  1
10  X<=Y?
11  GTO 00
12  RCL 01
13  RCL 02
14  P-R
15  1
16  -
17  STO 03
18  X<>Y
19  D-R
20  STO 04 
21  LASTX
22  R-D
23  RCL 00         
24  +
25  RCL 01
26  -
27  GTO 02
28  LBL 00
29  RCL 00
30  RCL 01
31  +
32  RCL 02
33  /
34  STO 04 
35  ENTER^
36  X^2
37  1
38  +
39  ST/ 04
40  SQRT
41  STO 03         
42  +
43  LN
44  RCL 01
45  -
46  RCL 02
47  X^2
48  RCL 03
49  *
50  ST/ 04
51  LASTX
52  RCL 02 
53  *
54  1/X
55  1
56  -
57  STO 03         
58  R^
59  LBL 02
60  ST+ X
61  ST* 04
62  RCL 03
63  ENTER^
64  X^2
65  RCL 04         
66  +
67  ABS
68  SQRT
69  -
70  /
71  ST- 01
72  ABS
73   E-6
74  X<Y?
75  GTO 01
76  RCL 01
77  END

 
   ( 95 bytes / SIZE 005 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             e            /
           X            M            E

 
Examples:

  e = 0.999   ENTER^
 M = 7°       XEQ "KPL"  >>>>  E = 52°27026151         ---Execution time = 10s---

  e = 1.2     ENTER^
 M = 0.1    XEQ "KPL"   >>>>  E = 0.423409616          ---Execution time = 7s---
 

     b) M-Code Routine
 

-This routine does not check for alpha data
 

08C  "L"
010   "P"
00B  "K"
2A0   SETDEC
0F8   READ 3(X)
128   WRIT 4(L)
2EE   ?C#0 ALL
3A0   ?NC RTN
0B8   READ 2(Y)
00E    A=0 ALL               A
1BE   A=A-1 MS             =
35C   PT=12
162    A=A+1 @PT         -1
01D   ?NCXQ                 C=
060    1807                    A+C
104    CLRF 8
2FE   ?C#0 MS
163    JNC +2Ch=+44d
108    SETF 8
0B8   READ 2(Y)
10E   A=C ALL
0F8   READ 3(X)
0A8   WRIT 2(Y)
0AE   A<>C ALL
0E8   WRIT 3(X)
070    N=C ALL           performs
3C4   ST=0                      the
1D5   ?NCXQ         polar-rectangular          ( in degrees without changing the current angular mode )   The results are in CPU registers C & N
078    1E75                  conversion
00E    A=0 ALL               A
1BE   A=A-1 MS             =
35C   PT=12
162    A=A+1 @PT         -1
01D   ?NCXQ                 C=
060    1807                    A+C
068    WRIT 1(Z)
0B0   C=N ALL
268   WRIT 9(Q)
205   ?NCXQ             performs D-R
064   1981                  provided N=C
10E   A=C ALL
278   READ 9(Q)
0AE  A<>C ALL
028   WRIT 0(T)
0AE  A<>C ALL
231   ?NCXQ             performs the
064   198C              R-D  conversion
0EE   C<>B ALL
0B8   READ 2(Y)
10E   A=C ALL
0F8   READ 3(X)
0A8   WRIT 2(Y)
0AE   A<>C ALL
0E8   WRIT 3(X)
2BE   C=-C-1 MS        C=-C
06E   A<>B ALL
01D  ?NCXQ                 C=
060   1807                    A+C
10E   A=C ALL
138   READ 4(L)
01B   JNC +03
1F3   JNC +3Eh=+62d
2AB  JNC -2Bh=-43d
01D  ?NCXQ                 C=
060   1807                    A+C
10E   A=C ALL
01D  ?NCXQ                 C=
060   1807                    A+C
070   N=C ALL
10E   A=C ALL
046   C=0 S&X               C
270   RAMSLCT             =
038   READDATA          T
135   ?NCXQ                 C=
060   184D                    A*C
268   WRIT 9(Q)
078   READ 1(Z)
10E   A=C ALL
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
278   READ 9(Q)
01D  ?NCXQ                 C=
060   1807                    A+C
05E   C=0 MS            C = | C |
2F9   ?NCXQ                C=
060    18BE                 sqrt(C)
2BE   C=-C-1 MS        C=-C
10E   A=C ALL
078   READ 1(Z)
01D  ?NCXQ                 C=
060   1807                    A+C
10E   A=C ALL
0B0   C=N ALL
2BE   C=-C-1 MS        C=-C
0AE   A<>C ALL
261    ?NCXQ               C=
060    1898                   A/C
070    N=C ALL
10E   A=C ALL
0F8   READ 3(X)
01D  ?NCXQ                 C=
060   1807                    A+C
0E8   WRIT 3(X)
0B0   C=N ALL
05E   C=0 MS            C = | C |
10E   A=C ALL
04E   C=0 ALL               C
2BE   C=-C-1 MS
35C   PT= 12
050    LD@PT- 1            =
21C   PT=2
250    LD@PT- 9
250    LD@PT- 9
110    LD@PT- 4       - 10^(-6)
01D  ?NCXQ                 C=
060   1807                    A+C
2FE   ?C#0 MS
360    ?C RTN
3CC   ?KEY                   press any key to stop an infinite loop
360    ?C RTN                or if the program is too slow
10C   ?FSET 8
227    JC -3Ch = -60d
0F8    READ 3(X)
10E    A=C ALL
013    JNC +02
20B   JNC -3Fh = -63d
138    READ 4(L)
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
0B8    READ 2(Y)
261    ?NCXQ               C=
060    1898                   A/C
070    N=C ALL
10E    A=C ALL
135   ?NCXQ                 C=
060   184D                    A*C
00E    A=0 ALL               A
35C   PT=12                    =
162    A=A+1 @PT         1
01D   ?NCXQ                 C=
060    1807                    A+C
068    WRIT 1(Z)
2F9   ?NCXQ                C=
060    18BE                 sqrt(C)
268    WRIT 9(Q)
10E    A=C ALL
078    READ 1(Z)
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
0B8   READ 2(Y)
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
0B8   READ 2(Y)
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
0B0   C=N ALL
0AE   A<>C ALL
261    ?NCXQ               C=
060    1898                   A/C
028    WRIT 0(T)
278    READ 9(Q)
10E    A=C ALL
0B8    READ 2(Y)
135    ?NCXQ                 C=
060    184D                    A*C
22D   ?NCXQ                C=
060    188B                    1/C
00E    A=0 ALL               A
1BE   A=A-1 MS             =
35C   PT=12
162    A=A+1 @PT         -1
01D   ?NCXQ                 C=
060    1807                    A+C
068    WRIT 1(Z)
278    READ 9(Q)
10E    A=C ALL
0B0   C=N ALL
01D   ?NCXQ                 C=
060    1807                    A+C
084    CLRF 5                 C
115    ?NCXQ                 =
06C   1B45                   Ln(C)
10E    A=C ALL
0F8    READ 3(X)
2BE   C=-C-1 MS        C=-C
207    JC -40h = -64d
 
 

      STACK        INPUTS      OUTPUTS
           Y             e            /
           X            M            E

 
Examples:

  e = 0.999   ENTER^
 M = 7°       XEQ "KPL"  >>>>  E = 52°27026151         ---Execution time = 5.5s---

  e = 1.2     ENTER^
 M = 0.1    XEQ "KPL"   >>>>  E = 0.423409616          ---Execution time = 3.1s---
 

Note:

-Register T = - f ''
-Register Z =  f '
 

4°)  Approximate Formula  ( e < 1 )
 

-Jean Meeus gives the formula:  tan E = ( sin M ) / ( cos M - e )
-The result is only an approximation but is still very good for the Sun ( maximum error = 0.2 arcsecond )
-For Pluto however ( e = 0.256 ), the error can reach  0.16°
 
 

 01  LBL "KPL"
 02  DEG
 03  1
 04  P-R
 05  RCL Z
 06  -
 07  R-P
 08  RDN
 09  END

 
( 18 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             e            e
           X            M            E

 
Example:

  e = 0.25   ENTER^
 M = 82°    XEQ "KPL"  >>>>  E = 96°3857  ( exact result = 96°2391045 )
 

References:

[1]   Jean Meeus - "Astronomical Algorithms" - Willmann-Bell  -  ISBN 0-943396-35-2
[2]   Robin M. Green - "Spherical Astronomy" - Cambridge University Press - ISBN  0-521-31779-7
[3]   http://mathworld.wolfram.com/KeplersEquation.html