Chebyshev-Eph

# 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  & 

-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 
-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
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
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
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
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
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
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
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
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:

   Newhall - "Numerical Representation of Planetary Ephemerides"
   ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii
   http://www.imcce.fr/inpop
   Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4