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