Ephemerides & Chebyshev Polynomials for the HP-41
Overview
1°) Positions
a) Recurrence Relation
b) Trigonometric Relation
c) Clenshaw's Algorithm
2°) Positions & Velocities
a) Recurrence Relation
b) Trigonometric Relation
c) Clenshaw's Algorithm
3°) M-Code Routines
a) Positions
b) Positions & Velocities
c) Positions & Velocities ( Clenshaw's
Algorithm )
d) Clenshaw's Algorithm ( simplified
version )
-One program is already listed in "Orthogonal Polynomials for the HP-41"
paragraph 5-b)
-Several variants are proposed hereunder.
-If y(t) is a coordinate of a celestial body at a given instant t ( t belonging to the interval [ t0 ; t0 + DT ] ) , y is given by the formula:
y = a0 + a1.T1(x) + ....... + an.Tn(x) where x = -1 + 2(t-t0)/DT ( -1 <= x <= +1 ) , Tk(x) are Chebyshev polynomials of the 1st kind
and a0 , a1 , ....... , an are given in various publications, see for example references [2] & [3]
-The derivative of such a function may also be obtained by
dy/dx = a1.U0(x) + 2a2.U1(x)
+........ + n.an.Un-1(x)
where Uk(x) are Chebyshev polynomials of the 2nd
kind
1°) Positions
a) Recurrence Relation
-Here, the Chebyshev polynomials of the 1st kind are computed by the
formula
Tk(x) = 2x.Tk-1(x)
- Tk-2(x) ; T0(x) = 1 ; T1(x)
= x
Data Registers: ( Registers Rbb thru Ree are to be initialized before executing "CdT" )
• Rbb = a0 • Rbb+1 = a1
......................... • Ree = an
ee = bb + n
Flags: /
Subroutines: /
| 01 LBL "CdT"
02 X<>Y 03 STO M 04 SIGN 05 STO O 06 X<>Y 07 STO N |
08 RCL IND M
09 ISG M 10 RCL IND M 11 RCL Z 12 * 13 + 14 STO P |
15 CLX
16 ISG M 17 LBL 01 18 R^ 19 ST+ X 20 RCL N 21 ST* Y |
22 X<> O
23 - 24 STO N 25 RCL IND M 26 * 27 + 28 ISG M |
29 GTO 01
30 RCL P 31 + 32 CLA 33 END |
( 59 bytes )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | bbb.eee | x |
| X | x | y(x) |
where bbb.eee is the control number of coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
Example: Calculate the mean ecliptic longitude of Jupiter ( with respect to the standard ecliptic J2000.0 ) on 2004/07/07 16h41m TT
-In the French Ephemeris "La Connaissance des Temps", we find 6 coefficients valid between 2004/01/00 0h TT to 2004/12/33 0h TT ( DT = 368 days )
-Store these coefficients into, say R10 thru R15 ( control number = 10.015 )
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
-The number days since 2004/01/00 0h TT is 189.695138889 whence x = -1 + 189.695138889 / 184 = 0.030951842
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time = 2s---
Notes:
-Unlike the version of "CdT" listed in "Orthogonal Polynomials for the
HP-41" paragraph 5-b)
this routine adds a0 + a1.x after
calculating a2.T2(x) + ....... + an.Tn(x)
-Since the first 2 coefficients usually have larger magnitudes than
the other ones,
we can expect less roundoff-errors in the final result.
-For the rest, the methods are identical.
b) Trigonometric Relation
-Here, we use Tk(x) = Cos ( k Acos x )
Data Registers: ( Registers R00 thru Rnn are to be initialized before executing "CdT2" )
• R00 = a0 • R01 = a1
......................... • Rnn = an
Flags: /
Subroutines: /
| 01 LBL "CdT2"
02 X<>Y 03 STO M 04 CLX 05 LBL 01 06 RCL Y 07 ACOS 08 RCL M 09 * 10 COS 11 RCL IND M 12 * 13 + 14 DSE M 15 GTO 01 16 RCL 00 17 + 18 END |
( 33 bytes / SIZE nnn+1 )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | n | x |
| X | x | y(x) |
where (n+1) is the number of coefficients a0
, a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
Example: The same one. Store the 6 coefficients into R00 thru R05
R00 = a0 = 173.010953
R01 = a1 = 13.996747
R02 = a2 = -0.032139
R03 = a3 = 0.003368
R04 = a4 = 0.000037
R05 = a5 = -0.000008
5
ENTER^
0.030951842 XEQ "CdT2"
>>>> L = 173°475979
---Execution time = 7s---
Notes:
-This version is obviously slower, but the routine is also shorter.
-Synthetic register M is cleared.
c) Clenshaw's Algorithm
-Another method to compute a sum of Chebyshev polynomials is given in
reference [4]
-Let:
• dn+1 = dn = 0
• dj = 2.x dj+1 - dj+2
+ aj for
j = n , n-1 , n-2 , ..... , 1
-Then y = a0 + a1.T1(x)
+ ....... + an.Tn(x) = x d1 - d2
+ a0
Data Registers: ( Registers Rbb thru Ree are to be initialized before executing "CdT" )
• Rbb = a0 • Rbb+1 = a1
......................... • Ree = an
ee = bb + n
Flags: /
Subroutines: /
| 01 LBL "CdT"
02 STO M 03 X<>Y 04 ENTER^ 05 INT 06 ST- Y 07 E3 08 ST* Z 09 / 10 + 11 STO N 12 CLST 13 LBL 01 14 RCL X 15 RCL M 16 ST+ X 17 * 18 R^ 19 - 20 RCL IND N 21 + 22 DSE N 23 GTO 01 24 RCL M 25 * 26 X<>Y 27 - 28 RCL IND N 29 + 30 CLA 31 END |
( 52 bytes )
| STACK | INPUTS | OUTPUTS |
| Y | bbb.eee | / |
| X | x | y(x) |
where bbb.eee is the control number of coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
Example: The same one. Store the 6 coefficients, for instance into R10 thru R15 ( control number = 10.015 )
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time = 2s---
Notes:
-The method is faster than the elementary one.
-The alpha register is cleared.
-If you place in Y a control number of the form eee.bbb
( like 15.010 in the example above ) , delete lines 04 to 10.
2°) Positions & Velocities
a) Recurrence Relations
-The Chebyshev polynomials of the 1st & 2nd kind are computed by
the formulae:
Tk(x) = 2x.Tk-1(x) - Tk-2(x)
; T0(x) = 1 ; T1(x) = x
( first kind )
Uk(x) = 2x.Uk-1(x) - Uk-2(x)
; U0(x) = 1 ; U1(x) = 2x
( second kind )
Data Registers: ( Registers Rbb thru Ree are to be initialized before executing "CdT" )
• Rbb = a0 • Rbb+1 = a1 ......................... • Ree = an ee = bb + n
Flag: F01 CF 01 = y(x)
SF 01 = dy/dx
Subroutines: /
| 01 LBL "CdT"
02 CLA 03 STO O 04 X<>Y 05 STO M 06 SIGN 07 FC? 01 08 STO N |
09 CHS
10 FS? 01 11 STO O 12 CLX 13 RCL IND M 14 FS? 01 15 CLX 16 ISG M |
17 LBL 01
18 RCL Y 19 ST+ X 20 RCL N 21 ST* Y 22 X<> O 23 - 24 STO N |
25 RCL IND M
26 * 27 ISG P 28 CLX 29 RCL P 30 FC? 01 31 SIGN 32 * |
33 +
34 ISG M 35 GTO 01 36 CLA 37 END |
( 66 bytes )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | bbb.eee | x |
| X | x | F(x) |
where bbb.eee is the control number of the coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
--> F(x) = y(x) if CF 01
--> F(x) = dy/dx if SF 01
Example: The same one. If the 6 coefficients are stored in R10 thru R15 , Control number = 10.015
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
• Position CF 01
CF 01
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time = 3s---
• Velocity SF 01
SF 01
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> dL/dx = 13.982645
---Execution time = 3s---
Note:
-To get the velocity in degree/day, the last result must be divided
by DT/2 , here 368/2 = 184
-It yields: dL/dt = 0.075992635 deg/day
b) Trigonometric Relation
-The Chebyshev polynomials of the 2nd kind are computed by
Uk(x) = Sin ( (k+1) Acos x ) / Sin ( Acos
x )
Data Registers: ( Registers R00 thru Rnn are to be initialized before executing "CdT2" )
• R00 = a0 • R01 = a1 ......................... • Rnn = an
Flag: F01
CF 01 = y(x)
SF 01 = dy/dx
Subroutines: /
| 01 LBL "CdT2"
02 DEG 03 X<>Y 04 STO M 05 CLX 06 LBL 01 07 RCL Y |
08 ACOS
09 RCL M 10 * 11 FS? 01 12 SIN 13 FC? 01 14 COS |
15 RCL M
16 FC? 01 17 SIGN 18 * 19 RCL IND M 20 * 21 + |
22 DSE M
23 GTO 01 24 X<>Y 25 ACOS 26 SIN 27 FC? 01 28 SIGN |
29 /
30 RCL 00 31 FS? 01 32 CLX 33 + 34 END |
( 55 bytes / SIZE nnn+1 )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | n | x |
| X | x | F(x) |
where (n+1) is the number of coefficients a0
, a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
--> F(x) = y(x) if CF 01
--> F(x) = dy/dx if SF 01
Example: The same one. Store the 6 coefficients into R00 thru R05
R00 = a0 = 173.010953
R01 = a1 = 13.996747
R02 = a2 = -0.032139
R03 = a3 = 0.003368
R04 = a4 = 0.000037
R05 = a5 = -0.000008
• Position CF 01
CF 01
5
ENTER^
0.030951842 XEQ "CdT2"
>>>> L = 173°475979
---Execution time = 9s---
• Velocity SF 01
SF 01
5
ENTER^
0.030951842 XEQ "CdT2"
>>>> dL/dx = 13.982645
---Execution time = 9s---
Notes:
-To get the velocity in degree/day, the last result must be divided
by DT/2 , here 368/2 = 184
-It yields: dL/dt = 0.075992635 deg/day
-This variant has a defect: if flag F01 is set, it does not work for
x = +/-1
c) Clenshaw's Algorithm
y(x) is obtained by the method of paragraph 1°-c)
To get the first derivative, let:
• d'n+1 = d'n = 0
• d'j = 2.x d'j+1 - d'j+2
+ j.aj for
j = n , n-1 , n-2 , ..... , 1
-Then d'1 = dy/dx = a1.U0(x)
+ 2a2.U1(x) +........ + n.an.Un-1(x)
Data Registers: ( Registers Rbb thru Ree are to be initialized before executing "CdT" )
• Rbb = a0 • Rbb+1 = a1 ......................... • Ree = an ee = bb + n
Flag: F01 CF 01 = y(x)
SF 01 = dy/dx
Subroutines: /
-Line 28 = TEXT0 or another NOP instruction like STO X
or LBL 00 ....
| 01 LBL "CdT"
02 STO M 03 X<>Y 04 INT 05 ST- L 06 LASTX 07 E3 08 * 09 STO N |
10 X<>Y
11 - 12 STO O 13 CLST 14 LBL 01 15 RCL X 16 RCL M 17 ST+ X 18 * |
19 R^
20 - 21 RCL IND N 22 RCL O 23 FC? 01 24 SIGN 25 * 26 + 27 DSE N |
28 ""
29 DSE O 30 GTO 01 31 FS? 01 32 GTO 02 33 RCL M 34 * 35 X<>Y 36 - |
37 RCL IND N
38 + 39 LBL 02 40 CLA 41 END |
( 67 bytes )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | bbb.eee | x |
| X | x | F(x) |
where bbb.eee is the control number of the coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
--> F(x) = y(x) if CF 01
--> F(x) = dy/dx if SF 01
Example: The same one. If the 6 coefficients are stored in R10 thru R15 , Control number = 10.015
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
• Position CF 01
CF 01
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time = 3s---
• Velocity SF 01
SF 01
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> dL/dx = 13.982645
---Execution time = 3s---
Note:
-To get the velocity in degree/day, the last result must be divided
by DT/2 , here 368/2 = 184
-It yields: dL/dt = 0.075992635 deg/day
3°) M-Code Routines
a) Positions
-This routine employs Clenshaw's algorithm ( paragraph 1-c) )
-13-digit routines are used to get a better precision.
-The existence of the registers is tested, but there is no check for
alpha data.
094 "T"
044 "d"
003 "C"
2A0 SETDEC
0B8 C=Y
070 N=C ALL
084 C
0ED =
064 frc(C)
226 C
226 =
226 1000.C
0F0 C<>N ALL
260 SETHEX
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
1BC RCR 11
0F0 C<>N ALL
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
106 A=C S&X
0B0 C=N ALL
11A A=C M
378 C=c
03C RCR 3
0E6 B<>C S&X
378 C=c
0E6 B<>C S&X
20E C=A+C ALL
23A C=C+1 M
070 N=C ALL
CPU register N now contains 1+bbb.eee ( absolute address in
hexadecimal )
106 A=C S&X
130 LDI S&X
200 200h
( correct value for an HP-41CV/CX or C with a QUAD memory module )
306 ?A<C S&X
381 ?NCGO
00A 02E0
NONEXISTENT is displayed if the specified register doesn't exist and the
routine stops
0F8 C=X
128 L=C
345 ?NCXQ
040 CLA
2A0 SETDEC
LOOP
138 C=L
10E A=C ALL
01D C=
060 A+C
1B8 C=N
158 M=C ALL
178 C=M
149 C=
060 AB*CM
238 C=P
2BE C=-C
158 M=C ALL
1F8 C=O
031 C=
060 AB+CM
178 C=M
1E8 O=C
1B8 C=N
228 P=C
0B0 C=N ALL
270 RAMSLCT
038 READATA
025 C=
060 AB+C
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
1A8 N=C
0CE C=B ALL
168 M=C
260 SETHEX
0B0 C=N ALL
266 C=C-1 S&X
070 N=C ALL
106 A=C S&X
03C RCR 3
306 ?A<C S&X
2D3 JNC-38d
GOTO LOOP
2A0 SETDEC
1B8 C=N
10E A=C ALL
138 C=L
13D C=
060 AB*C
238 C=P
2BE C=-C
158 M=C ALL
1F8 C=O
031 C=
060 AB+CM
0B0 C=N ALL
270 RAMSLCT
038 READATA
025 C=
060 AB+C
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
0E8 X=C
345 ?NCGO
042 CLA
( 104 words )
| STACK | INPUTS | OUTPUTS |
| T | T | T |
| Z | Z | Z |
| Y | bbb.eee | bbb.eee |
| X | x | y(x) |
| L | / | x |
where bbb.eee is the control number of the coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
Example: Always the same one. If the 6 coefficients are stored in R10 thru R15 , Control number = 10.015
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time < 1s---
Notes:
-The alpha register is cleared but synthetic register Q is not
used.
-If eee = 000 , you'll get a DATA ERROR message, even in R00
does exist.
-This routine is the fastest program that is written in this page.
-With 14 coefficients, execution time = 1.4 second
-With 100 coefficients, ------------- = 9.6 seconds.
b) Positions & Velocities
-This M-Code routine "CdT" employs the recurrence relations of paragraphs
1-a) & 2-a)
-13-digit routines are used to get a better precision.
-The existence of the registers is tested, but there is no check for
alpha data.
094 "T"
044 "d"
003 "C"
2A0 SETDEC
0B8 C=Y
070 N=C ALL
084 C
0ED =
064 frc(C)
226 C
226 =
226 1000.C
0F0 C<>N ALL
260 SETHEX
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
1BC RCR 11
0F0 C<>N ALL
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
106 A=C S&X
0B0 C=N ALL
11A A=C M
378 C=c
03C RCR 3
0E6 B<>C S&X
378 C=c
0E6 B<>C S&X
20E C=A+C ALL
070 N=C ALL
CPU register N now contains bbb.eee ( absolute address in
hexadecimal )
106 A=C S&X
03C RCR 3
306 ?A<C S&X
013 JNC+02
0A6 A<>C S&X
130 LDI S&X
200 200h
( correct value for an HP-41CV/CX or C with a QUAD memory module )
306 ?A<C S&X
381 ?NCGO
00A 02E0
NONEXISTENT is displayed if the specified register doesn't exist and the
routine stops
2A0 SETDEC
0F8 C=X
128 L=C
104 CLRF 8
0B8 C=Y
2FE ?C<0
013 JNC+02
108 SETF 8
04E C=0 ALL
0E8 X=C
1A8 N=C
10C ?FSET 8
053 JNC+10d
168 M=C
35C C=
050 1
1E8 O=C
04E C=0 ALL
2BE C=-C
228 P=C
1A0 A=B=C=0
08B JNC+17d
35C C=
050 1
168 M=C
138 C=L
02E B=0 ALL
0FA B<>C M
228 P=C
0EE B<>C ALL
1E8 O=C
0B0 C=N ALL
03C RCR 3
270 RAMSLCT
038 READATA
02E C
0FA ->
0AE AB
081 0 ramslct
064 AB STO Q+
260 SETHEX
LOOP
0B0 C=N ALL
23A C=C+1 M
070 N=C ALL
106 A=C S&X
03C RCR 3
306 ?A<C S&X
1BF JC+55d
2A0 SETDEC
138 C=L
10E A=C ALL
01D C=
060 A+C
1B8 C=N
158 M=C ALL
178 C=M
149 C=
060 AB*CM
238 C=P
2BE C=-C
158 M=C ALL
1F8 C=O
031 C=
060 AB+CM
178 C=M
1E8 O=C
1B8 C=N
228 P=C
0AE A<>C ALL
1A8 N=C
0EE B<>C ALL
168 M=C
0F8 C=X
10C ?FSET 8
017 JC+02
1A0 A=B=C=0
02E C
0FA ->
0AE AB
001 C=
060 AB+1
0E8 X=C
0B0 C=N ALL
03C RCR 3
270 RAMSLCT
038 READATA
13D C=
060 AB*C
046 C=0 S&X
270 RAMSLCT
1B8 C=N
158 M=C ALL
178 C=M
149 C=
060 AB*CM
0D1 RCL
064 Q+
031 C=
060 AB+CM
089 AB
064 STO Q+
21B JNC-61d
GOTO LOOP
2A0 SETDEC
1A0 A=B=C=0
0D1 RCL
064 Q+
031 C=
060 AB+CM
0E8 X=C
345 ?NCGO
042 CLA
( 151 words )
| STACK | INPUTS | OUTPUTS |
| T | T | T |
| Z | Z | Z |
| Y | +/- bbb.eee | +/- bbb.eee |
| X | x | F(x) |
| L | / | x |
where bbb.eee is the control number of the coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
--> F(x) = y(x) if the control number is
+ bbb.eee
--> F(x) = dy/dx if the control number is -
bbb.eee
Example: Always the same one. If the 6 coefficients are stored in R10 thru R15 , Control number = 10.015
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
• Position:
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time = 1s---
• Velocity:
10.015 CHS ENTER^
0.030951842 XEQ "CdT"
>>>> dL/dx = 13.982645
---Execution time = 1s---
Notes:
-To get the velocity in degree/day, the last result must be divided
by DT/2 , here 368/2 = 184
-It yields: dL/dt = 0.075992635 deg/day
-The alpha register is cleared and synthetic register Q is used.
-If eee = 000 , you'll get a DATA ERROR message, even if R00
does exist.
c) Positions & Velocities
( Clenshaw's Algorithm )
-This routine employs Clenshaw's algorithm ( paragraphs 1-c) & 2-c)
)
-13-digit routines are used to get a better precision.
-The existence of the registers is tested, but there is no check for
alpha data.
094 "T"
044 "d"
003 "C"
2A0 SETDEC
0B8 C=Y
070 N=C ALL
084 C
0ED =
064 frc(C)
226 C
226 =
226 1000.C
268 Q=C
0F0 C<>N ALL
260 SETHEX
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
1BC RCR 11
0F0 C<>N ALL
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
106 A=C S&X
0B0 C=N ALL
11A A=C M
378 C=c
03C RCR 3
0E6 B<>C S&X
378 C=c
0E6 B<>C S&X
20E C=A+C ALL
23A C=C+1 M
070 N=C ALL
CPU register N now contains 1+bbb.eee ( absolute address in
hexadecimal )
106 A=C S&X
130 LDI S&X
200 200h
( correct value for an HP-41CV/CX or C with a QUAD memory module )
306 ?A<C S&X
381 ?NCGO
00A 02E0
NONEXISTENT is displayed if the specified register doesn't exist and the
routine stops
0F8 C=X
128 L=C
345 ?NCXQ
040 CLA
2A0 SETDEC
278 C=Q
088 C
0ED =
064 Int(C)
268 Q=C
0B8 C=Y
088 C
0ED =
064 Int(C)
2BE C=-C
10E A=C ALL
278 C=Q
01D C=
060 A+C
05E C= | C |
268 Q=C
104 CLRF 8
0B8 C=Y
2FE ?C<0
013 JNC+02
108 SETF 8
2A0 SETDEC
LOOP
138 C=L
10E A=C ALL
01D C=
060 A+C
1B8 C=N
158 M=C ALL
178 C=M
149 C=
060 AB*CM
238 C=P
2BE C=-C
158 M=C ALL
1F8 C=O
031 C=
060 AB+CM
178 C=M
1E8 O=C
1B8 C=N
228 P=C
0AE A<>C ALL
1A8 N=C
0CE C=B ALL
168 M=C
278 C=Q
10C ?FSET 8
027 JC+04
04E C
35C =
050 1
10E A=C ALL
0B0 C=N ALL
270 RAMSLCT
038 READATA
135 C=
060 A*C
046 C=0 S&X
270 RAMSLCT
1B8 C=N
158 M=C ALL
178 C=M
031 C=
060 AB+CM
0AE A<>C ALL
1A8 N=C
0CE C=B ALL
168 M=C
278 C=Q
02E C
0FA ->
0AE AB
009 C=
060 AB-1
268 Q=C
260 SETHEX
0B0 C=N ALL
266 C=C-1 S&X
070 N=C ALL
106 A=C S&X
03C RCR 3
306 ?A<C S&X
21B JNC-61d
GOTO LOOP
2A0 SETDEC
1B8 C=N
10E A=C ALL
178 C=M
0EE B<>C ALL
04E C=0 ALL
10C ?FSET 8
06F JC+13d
138 C=L
13D C=
060 AB*C
238 C=P
2BE C=-C
158 M=C ALL
1F8 C=O
031 C=
060 AB+CM
0B0 C=N ALL
270 RAMSLCT
038 READATA
025 C=
060 AB+C
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
0E8 X=C
345 ?NCGO
042 CLA
( 155 words )
| STACK | INPUTS | OUTPUTS |
| T | T | T |
| Z | Z | Z |
| Y | +/- bbb.eee | +/- bbb.eee |
| X | x | F(x) |
| L | / | x |
where bbb.eee is the control number of the coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
--> F(x) = y(x) if the control number is
+ bbb.eee
--> F(x) = dy/dx if the control number is -
bbb.eee
Example: Always the same one. If the 6 coefficients are stored in R10 thru R15 , Control number = 10.015
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
• Position:
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time = 1s---
• Velocity:
10.015 CHS ENTER^
0.030951842 XEQ "CdT"
>>>> dL/dx = 13.982645
---Execution time = 1s---
Notes:
-To get the velocity in degree/day, the last result must be divided
by DT/2 , here 368/2 = 184
-It yields: dL/dt = 0.075992635 deg/day
-The alpha register is cleared and synthetic register Q is used.
-If eee = 000 , you'll get a DATA ERROR message, even in R00
does exist.
-With 14 coefficients, execution time = 2.2 seconds
-With 100 coefficients, ------------- = 16 seconds.
-All these M-Code routines use 2 registers for the dj's:
1 for the 13-digit mantissa and 1 for the sign & exponent.
-So, they could be simplified with perhaps a small roundoff-error in
the last decimal...
d) Clenshaw's Algorithm
( simplified version )
-This M-code routine uses only one register for each intermediate result:
-The program runs faster with perhaps a minor loss of precision.
-The existence of the registers is tested, but there is no check for
alpha data.
094 "T"
044 "d"
003 "C"
2A0 SETDEC
0B8 C=Y
070 N=C ALL
084 C
0ED =
064 frc(C)
226 C
226 =
226 1000.C
268 Q=C
0F0 C<>N ALL
260 SETHEX
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
1BC RCR 11
0F0 C<>N ALL
38D ?NCXQ
converts the decimal number in C
008 02E3
into hexadecimal in C S&X field
106 A=C S&X
0B0 C=N ALL
11A A=C M
378 C=c
15A A=A+C M
03C RCR 3
146 A=A+C S&X
0AE A<>C ALL
23A C=C+1 M
070 N=C ALL
CPU register N now contains 1+bbb.eee ( absolute address in
hexadecimal )
106 A=C S&X
130 LDI S&X
200 200h
( correct value for an HP-41CV/CX or C with a QUAD memory module )
306 ?A<C S&X
381 ?NCGO
00A 02E0
NONEXISTENT is displayed if the specified register doesn't exist and the
routine stops
0F8 C=X
128 L=C
345 ?NCXQ
040 CLA
2A0 SETDEC
278 C=Q
088 C
0ED =
064 Int(C)
268 Q=C
0B8 C=Y
088 C
0ED =
064 Int(C)
2BE C=-C
10E A=C ALL
278 C=Q
01D C=
060 A+C
05E C= | C |
268 Q=C
104 CLRF 8
0B8 C=Y
2FE ?C<0
013 JNC+02
108 SETF 8
2A0 SETDEC
LOOP
138 C=L
10E A=C ALL
01D C=
060 A+C
178 C=M
13D C=
060 AB*C
1B8 C=N
2BE C=-C
025 C=
061 AB+C
1E8 O=C
178 C=M
1A8 N=C
0B0 C=N ALL
270 RAMSLCT
038 READATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
10C ?FSET 8
02B JNC+05
278 C=Q
135 C=
060 A*C
10E A=C ALL
1F8 C=O
01D C=
060 A+C
168 M=C
10C ?FSET 8
043 JNC+08
278 C=Q
02E C
0FA ->
0AE AB
009 C=
060 AB-1
268 Q=C
260 SETHEX
0B0 C=N ALL
266 C=C-1 S&X
070 N=C ALL
106 A=C S&X
03C RCR 3
306 ?A<C S&X
28B JNC-47d
GOTO LOOP
2A0 SETDEC
178 C=M
10C ?FSET 8
097 JC+18d
10E A=C ALL
138 C=L
135 C=
060 A*C
1B8 C=N
2BE C=-C
025 C=
061 AB+C
0B0 C=N ALL
270 RAMSLCT
038 READATA
025 C=
060 AB+C
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
0E8 X=C
345 ?NCGO
042 CLA
( 135 words )
| STACK | INPUTS | OUTPUTS |
| T | T | T |
| Z | Z | Z |
| Y | +/- bbb.eee | +/- bbb.eee |
| X | x | F(x) |
| L | / | x |
where bbb.eee is the control number of the coefficients
a0 , a1 , ....... , an
and x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
--> F(x) = y(x) if the control number is
+ bbb.eee
--> F(x) = dy/dx if the control number is -
bbb.eee
Example: Again the same one. If the 6 coefficients are stored in R10 thru R15 , Control number = 10.015
R10 = a0 = 173.010953
R11 = a1 = 13.996747
R12 = a2 = -0.032139
R13 = a3 = 0.003368
R14 = a4 = 0.000037
R15 = a5 = -0.000008
• Position:
10.015
ENTER^
0.030951842 XEQ "CdT"
>>>> L = 173°475979
---Execution time < 1s---
• Velocity:
10.015 CHS ENTER^
0.030951842 XEQ "CdT"
>>>> dL/dx = 13.982645
---Execution time < 1s---
Notes:
-To get the velocity in degree/day, the last result must be divided
by DT/2 , here 368/2 = 184
-It yields: dL/dt = 0.075992635 deg/day
-The alpha register is cleared and synthetic register Q is used.
-If eee = 000 , you'll get a DATA ERROR message, even in R00
does exist.
-With 14 coefficients, execution time = 1.3 second
-With 100 coefficients, ------------- = 10 seconds ( 14s for
dy/dx ).
References:
[1] Newhall - "Numerical Representation of Planetary Ephemerides"
[2] ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii
[3] http://www.imcce.fr/inpop
[4] Press , Vetterling , Teukolsky , Flannery - "Numerical
Recipes in C++" - Cambridge University Press - ISBN 0-521-75033-4