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
128   WRIT 4(L)
2EE   ?C#0 ALL
3A0   ?NC RTN
068   WRIT 1(Z)
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
231   ?NCXQ             performs the
064   198C              R-D  conversion
028   WRIT 0(T)
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          =
135   ?NCXQ              C=
060   184D                 A*C
10E   A=C ALL
01D  ?NCXQ              C=
060   1807                  A+C
0EE   C<>B ALL
10E   A=C ALL
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
10E   A=C ALL
01D  ?NCXQ              C=
060   1807                  A+C
10E   A=C ALL
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
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

( 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
10E   A=C ALL
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
10E    A=C ALL
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
10E    A=C ALL
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
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
261    ?NCXQ                C=
060    1898                    A/C
070    N=C ALL
10E    A=C ALL
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             =
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
10E    A=C ALL
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
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
2BE   C=-C-1 MS        C=-C
01D   ?CXQ                   C=
061    1807                    A+C
068    WRIT 1(Z)
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
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
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

( 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
128   WRIT 4(L)
2EE   ?C#0 ALL
3A0   ?NC RTN
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
10E   A=C ALL
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
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
10E   A=C ALL
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
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             =
135   ?NCXQ                 C=
060   184D                    A*C
268   WRIT 9(Q)
10E   A=C ALL
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
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
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
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
10E    A=C ALL
013    JNC +02
20B   JNC -3Fh = -63d
01D   ?NCXQ                 C=
060    1807                    A+C
10E    A=C ALL
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
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
135   ?NCXQ                 C=
060   184D                    A*C
10E   A=C ALL
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)
10E    A=C ALL
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)
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
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