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  [ t₀ ; t₀ + DT ] ) ,  y is given by the formula:

       y = a₀ + a₁.T₁(x) + ....... + a_n.T_n(x)             where  x = -1 + 2(t-t₀)/DT      ( -1 <= x <= +1 )  ,   T_k(x) are Chebyshev polynomials of the 1st kind

  and    a₀ , a₁ , ....... , a_n  are given in various publications, see for example references [2] & [3]

-The derivative of such a function may also be obtained by

     dy/dx = a₁.U₀(x) + 2a₂.U₁(x) +........ + n.a_n.U_n-1(x)      where   U_k(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
 

      T_k(x) = 2x.T_k-1(x) - T_k-2(x)  ;  T₀(x) = 1  ;  T₁(x) = x
 

Data Registers:               ( Registers Rbb thru Ree are to be initialized before executing "CdT" )

                                            •  Rbb = a₀    •  Rbb+1 = a₁    .........................    •  Ree = a_n        ee = bb + n
Flags: /
Subroutines: /
 
 

 
     ( 59 bytes )
 
 

  where  bbb.eee  is the control number of coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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 a₀ + a₁.x  after   calculating  a₂.T₂(x) + ....... + a_n.T_n(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   T_k(x) = Cos ( k Acos x )
 

Data Registers:                ( Registers R00 thru Rnn are to be initialized before executing "CdT2" )

                                            •  R00 = a₀    •  R01 = a₁    .........................    •  Rnn = a_n
Flags: /
Subroutines: /
 
 

 
( 33 bytes / SIZE nnn+1 )
 
 

  where  (n+1) is the number of coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/DT      ( -1 <= x <= +1 )

Example:    The same one. Store the 6 coefficients into R00 thru R05

    R00 = a₀ = 173.010953
    R01 = a₁ =   13.996747
    R02 = a₂ =    -0.032139
    R03 = a₃ =     0.003368
    R04 = a₄ =     0.000037
    R05 = a₅ =   -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:

  •  d_n+1 = d_n = 0
  •  d_j = 2.x d_j+1 - d_j+2 + a_j        for   j = n , n-1 , n-2 , ..... , 1

-Then   y = a₀ + a₁.T₁(x) + ....... + a_n.T_n(x)  = x d₁ - d₂ + a₀
 

Data Registers:                ( Registers Rbb thru Ree are to be initialized before executing "CdT" )

                                            •  Rbb = a₀    •  Rbb+1 = a₁    .........................    •  Ree = a_n        ee = bb + n
Flags: /
Subroutines: /
 
 

 
    ( 52 bytes )
 
 

  where  bbb.eee  is the control number of coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/DT      ( -1 <= x <= +1 )

Example:    The same one. Store the 6 coefficients, for instance into R10 thru R15  ( control number = 10.015 )

    R10 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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:
 

  T_k(x) = 2x.T_k-1(x) - T_k-2(x)  ;  T₀(x) = 1  ;  T₁(x) = x             ( first kind )
  U_k(x) = 2x.U_k-1(x) - U_k-2(x)  ;  U₀(x) = 1  ;  U₁(x) = 2x      ( second kind )
 

Data Registers:                ( Registers Rbb thru Ree are to be initialized before executing "CdT" )

                                            •  Rbb = a₀    •  Rbb+1 = a₁    .........................    •  Ree = a_n        ee = bb + n

Flag:  F01    CF 01 = y(x)
                     SF 01 = dy/dx
Subroutines: /
 
 

 
     ( 66 bytes )
 
 

  where  bbb.eee is the control number of the coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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
 

    U_k(x) = Sin ( (k+1) Acos x ) / Sin ( Acos x )
 

Data Registers:                ( Registers R00 thru Rnn are to be initialized before executing "CdT2" )

                                            •  R00 = a₀    •  R01 = a₁    .........................    •  Rnn = a_n

Flag:  F01          CF 01 = y(x)
                           SF 01 = dy/dx
Subroutines: /
 
 

 
    ( 55 bytes / SIZE nnn+1 )
 
 

  where  (n+1) is the number of coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R01 = a₁ =   13.996747
    R02 = a₂ =    -0.032139
    R03 = a₃ =     0.003368
    R04 = a₄ =     0.000037
    R05 = a₅ =   -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.a_j        for   j = n , n-1 , n-2 , ..... , 1

-Then   d'₁ = dy/dx = a₁.U₀(x) + 2a₂.U₁(x) +........ + n.a_n.U_n-1(x)
 

Data Registers:                ( Registers Rbb thru Ree are to be initialized before executing "CdT" )

                                            •  Rbb = a₀    •  Rbb+1 = a₁    .........................    •  Ree = a_n        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  ....
 
 

 
     ( 67 bytes )
 
 

  where  bbb.eee is the control number of the coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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 )
 
 

  where  bbb.eee is the control number of the coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/DT      ( -1 <= x <= +1 )

Example:    Always the same one. If the 6 coefficients are stored in R10 thru R15 ,   Control number = 10.015

    R10 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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 )
 
 

  where  bbb.eee is the control number of the coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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 )
 
 

  where  bbb.eee is the control number of the coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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 d_j'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 )
 
 

  where  bbb.eee is the control number of the coefficients   a₀ , a₁ , ....... , a_n
    and      x = -1 + 2(t-t₀)/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 = a₀ = 173.010953
    R11 = a₁ =   13.996747
    R12 = a₂ =    -0.032139
    R13 = a₃ =     0.003368
    R14 = a₄ =     0.000037
    R15 = a₅ =   -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