All code in this article is public domain

I have implemented an example implementation based on the description below to show how everything fits together.

The JPL Development Ephemeris provide high quality, accurate data on the positions of the planets, the moon, and sometimes a few other variables. Writing software to use the JPL files is not difficult. Unfortunately the documentation is pretty lacking, scattered, and mostly relies on reading the source code of reference implementations. Having looked at the reference implementations, I'm pretty convinced that the authors are better astronomers than they are programmers. I'm not here to judge, I'm certainly a much better programmer than I am an astronomer. I'm hoping with this article I can help bridge the gap, and explain how to use the JPL ephemeris data in terms that are easier for programmers to follow who might not be as informed on astronomy, and hopefully open up new opprotunities.

If you've looked at some of the implementations provided by NASA, or the Project Pluto implementations, it might look like it's an impossibly hard task to grapple. But the reality is that these are optimised implementations, designed to work with multiple ephemeris versions, and are written in pretty old languages that lack features that make code easier to understand (even if they run slower). In reality, in most of these implementations, only about 10-30 lines of code actualy deal with calculating the position (or velocity) of the planets. Most of the rest of the code is dealing with edge cases for dealing with multiple versions, memory management, and finding the right coefficients to finally pass to the computation function.

If you're willing to forego the flexibility of being able to process ephemeris files you've never seen before, then writing a program to compute all of the series in a JPL ephemeris version can be done in about 100 lines of code in pretty much any modern language. In fact, if you were to rearrange the JPL ephemeris files into a format which is easier to parse, you could likely get that down to 10 to 20 lines of code. The point I'm trying to make here, is that it's the file format that is complicated, not the actual computation of positions.

Probably the number one reason to write your own software, specific to your needs, is that you won't be locked in to using the format JPL provides the ephemeris data in. As you'll see below, much of the data is independent of each other. And, if you're only interested in data on the Sun, Earth, and Moon, you can extract that data into your own files, considerably reducing the size of the files you have to work with.

Another reason is most of the implementations recommended by the JPL are written in C, which is great if speed is your only concern. But modern computing environments aren't very friendly to the C ecosystem. A web browser won't let you run a C application at all, but is very likely the most widely used platform for applications. Platforms like Android will let you run C, but providing binaries for every type of device is practically impossible. Other systems, like microcontrollers, have no concept of file systems.

So, there's a lot of room to grow. So, let's get started.

KSIZE= 2036 NCOEFF= 1018 GROUP 1010 JPL Planetary Ephemeris DE405/DE405 Start Epoch: JED= 2305424.5 1599 DEC 09 00:00:00 Final Epoch: JED= 2525008.5 2201 FEB 20 00:00:00 GROUP 1030 2305424.50 2525008.50 32. GROUP 1040 156 DENUM LENUM TDATEF TDATEB CENTER CLIGHT AU EMRAT GM1 GM2 GMB GM4 GM5 GM6 GM7 GM8 GM9 GMS RAD1 RAD2 RAD4 JDEPOC X1 Y1 Z1 XD1 YD1 ZD1 X2 Y2 Z2 XD2 YD2 ZD2 XB YB ZB XDB YDB ZDB X4 Y4 Z4 XD4 YD4 ZD4 X5 Y5 Z5 XD5 YD5 ZD5 X6 Y6 Z6 XD6 YD6 ZD6 X7 Y7 Z7 XD7 YD7 ZD7 X8 Y8 Z8 XD8 YD8 ZD8 X9 Y9 Z9 XD9 YD9 ZD9 XM YM ZM XDM YDM ZDM XS YS ZS XDS YDS ZDS BETA GAMMA J2SUN GDOT MA0001 MA0002 MA0004 MAD1 MAD2 MAD3 RE ASUN PHI THT PSI OMEGAX OMEGAY OMEGAZ AM J2M J3M J4M C22M C31M C32M C33M S31M S32M S33M C41M C42M C43M C44M S41M S42M S43M S44M LBET LGAM K2M TAUM AE J2E J3E J4E K2E0 K2E1 K2E2 TAUE0 TAUE1 TAUE2 DROTEX DROTEY GMAST1 GMAST2 GMAST3 KVC IFAC PHIC THTC PSIC OMGCX OMGCY OMGCZ PSIDOT MGMIS ROTEX ROTEY GROUP 1041 156 0.405000000000000000D+03 0.405000000000000000D+03 0.000000000000000000D+00 0.119970525194723000D+17 0.000000000000000000D+00 0.299792457999999984D+06 0.149597870691000015D+09 0.813005600000000044D+02 0.491254745145081187D-10 0.724345248616270270D-09 0.899701134671249882D-09 0.954953510577925806D-10 0.282534590952422643D-06 0.845971518568065874D-07 0.129202491678196939D-07 0.152435890078427628D-07 0.218869976542596968D-11 0.295912208285591095D-03 0.243976000000000022D+04 0.605230000000000018D+04 0.339751499999999987D+04 0.244040050000000000D+07 0.361762714603509283D+00 -0.907819677295860494D-01 -0.857149831817633490D-01 0.336749391398414016D-02 0.248945204676488743D-01 .... 0.646682543384255465D-13 0.127748118910414607D-13 0.333405877296029502D-14 0.000000000000000000D+00 0.299999999999999974D-03 -0.425951830000000000D-02 0.408844299999999994D+00 -0.171450900000000006D+01 0.000000000000000000D+00 -0.158167070000000005D-05 0.229888000000000009D+00 0.000000000000000000D+00 0.100000000000000000D+01 0.000000000000000000D+00 0.000000000000000000D+00 GROUP 1050 3 171 231 309 342 366 387 405 423 441 753 819 899 14 10 13 11 8 7 6 6 6 13 11 10 10 4 2 2 1 1 1 1 1 1 8 2 4 4 GROUP 1070

The first two numbers here are Julian Days. They are the earliest, and latest dates that the ephemeris can generate data for. You'll notice that GROUP 1010, just above GROUP 1030 has the same data, and with their respective Gregorian Calendar dates listed along side. But if you're parsing the file automatically, you'll likely find GROUP 1030 easier to parse.

The last number, which is 32 in this case, is the total number of days a given block contains data for. As I said above, each file is broken up into many blocks, and each block is only good for a set number of days. You'll use this later to compute an offset into the full list of data to find the right coefficients.

This is the most complicated part of the data, and will be explained in more detail when we actually look at the data files. But the short version is that each column here represents a planet (or some other data element). The order is consistent across all of the ephemeris versions. The order is:

Table 1. List of Series order in the JPL Files:

# Properties | Units | Center | Name | |
---|---|---|---|---|

1 | 3 | km | SSB | Mercury |

2 | 3 | km | SSB | Venus |

3 | 3 | km | SSB | Earth-Moon barycenter |

4 | 3 | km | SSB | Mars |

5 | 3 | km | SSB | Jupiter |

6 | 3 | km | SSB | Saturn |

7 | 3 | km | SSB | Uranus |

8 | 3 | km | SSB | Neptune |

9 | 3 | km | SSB | Pluto |

10 | 3 | km | Earth | Moon (geocentric) |

11 | 3 | km | SSB | Sun |

12 | 2 | radians | Earth Nutations in longitude and obliquity (IAU 1980 model) | |

13 | 3 | radians | Lunar mantle libration | |

14 | 3 | radians/day | Lunar mantle angular velocity | |

15 | 1 | Seconds | TT-TDB (at geocenter) |

You'll notice that not all ephemeris versions contain data for every data element. For example DE405 (above) does not have data for Lunar Mantle Velocity, nor TT-TDB. and DE432 (below) does not have data for Nutations, Mantle Velocity, nor TT-TDB. Exactly how each ephemeris version expresses the fact that they lack data for a certain element is somewhat variable. You'll notice that DE405 just lacks those columns, where DE432 has indicated 0 for the number of coefficients. Otherwise, you'll note that the numbers are pretty much the same between DE405 and DE432.

GROUP 1050 3 171 231 309 342 366 387 405 423 441 753 819 819 939 939 14 10 13 11 8 7 6 6 6 13 11 0 10 0 0 4 2 2 1 1 1 1 1 1 8 2 0 4 0 0

The numbers in each of these columns (from top to bottom) are:

- The start offset of the series in a block
- The number of coefficients for each property
- The number of subintervals

Not included here, but also important is the number of properties for each series. This information you have to infer by the type of series. Though you could compute it by looking at where the next series starts (or for the last series, the length of a block). Each planet has three properties (X, Y, and Z). Nutations have two properties (dpsi, deps). Librations have three (phi, theta, psi). Mantel velocity has three (omega x, omega y, omega z). And TT-TDB has just one. This data is included in Table 1 above.

```
```**Earth = EMB - Moon / ( 1 + earthMoonMassRatio )**

So, anytime you need tht position of the Earth, you must first compute the EMB and Moon positions, then use those positions in the equation above. The earthMoonMassRatio value comes from the header file, and is highlighted above (0.813005600000000044D+02).

JPL provides ASCII versions of the ephemeris files, and binary versions of some of them, some in little-endian and some in big-endian format. For the purposes of this article I will deal only with the ASCII format. The binary format is pretty similar once you get past the header. As I said above, one of the biggest advantages to writing your own ephemeris generator is that you're not locked in to using the JPL formatted files. Converting them to binary can speed up processing significantly, and if that's important to your application, you can probably come up with a more appropriate binary format than the JPL version. For that reason, I won't go into the binary format any further.

The general format of the ASCII files, using DE405 as an example, is as follows

- Block
- Block Start Date - 1 Double
- Block End Date - 1 Double
- Mercury Series: Offset 3
- Mercury Subinterval 1
- X Coefficients - 14 Doubles
- Y Coefficients - 14 Doubles
- Z Coefficients - 14 Doubles
- Mercury Subinterval 2
- X Coefficients - 14 Doubles
- Y Coefficients - 14 Doubles
- Z Coefficients - 14 Doubles
- Mercury Subinterval 3
- X Coefficients - 14 Doubles
- Y Coefficients - 14 Doubles
- Z Coefficients - 14 Doubles
- Mercury Subinterval 4
- X Coefficients - 14 Doubles
- Y Coefficients - 14 Doubles
- Z Coefficients - 14 Doubles
- Venus Series: Offset 171
- Venus Subinterval 1
- X Coefficients - 10 Doubles
- Y Coefficients - 10 Doubles
- Z Coefficients - 10 Doubles
- Venus Subinterval 2
- X Coefficients - 10 Doubles
- Y Coefficients - 10 Doubles
- Z Coefficients - 10 Doubles
- Earth-Moon Barycenter Series: Offset 231
- EMB Subinterval 1
- X Coefficients - 13 Doubles
- Y Coefficients - 13 Doubles
- Z Coefficients - 13 Doubles
- EMB Subinterval 2
- X Coefficients - 13 Doubles
- Y Coefficients - 13 Doubles
- Z Coefficients - 13 Doubles
- . . . . . . .
- . . . . . . .
- And So On
- . . . . . . .
- . . . . . . .
- Lunar mantle libration: Offset 899
- Libration Subinterval 1
- Phi Φ Coefficients - 10 Doubles
- Theta θ Coefficients - 10 Doubles
- Psi ψ Coefficients - 10 Doubles
- Libration Subinterval 2
- Phi Φ Coefficients - 10 Doubles
- Theta θ Coefficients - 10 Doubles
- Psi ψ Coefficients - 10 Doubles
- Libration Subinterval 3
- Phi Φ Coefficients - 10 Doubles
- Theta θ Coefficients - 10 Doubles
- Psi ψ Coefficients - 10 Doubles
- Libration Subinterval 4
- Phi Φ Coefficients - 10 Doubles
- Theta θ Coefficients - 10 Doubles
- Psi ψ Coefficients - 10 Doubles

Each ephemeris version has files named slightly differently, mostly based on how many years of data each file contains. The basic format of the filename is:

```
ASC[P/M][Year].[Version]
```

Where [P/M] stand for "plus" or "minus". It is "M" for negative years in the Gregorian Calendar, and "P" otherwise. [Year] is the first year the
file has data for. And [Version] is the version of the ephemeris it was generated for. Since the 32 day blocks don't line up exactly on Gergorian
Calendar years, the start and end dates will have some overlap with other files. The first and last blocks of data for 32 days is repeated in
adjacent files. So, once you've computed the Gregorian Year you're looking for, you know which file to look in. But you need to be careful when
joining two files, otherwise including duplicated data will likey disrupt the math used to compute indexes into the data. The year is **normally** four
digits, but after DE430 some have 5 digits for the year.

The numbers in the files are formatted as exponentials, with a D denoting the exponent. You will need to replace this D with an E in order for it to parse correctly in most programming languages.

If you open up the ascp2020.405 file as an example. The first thing you'll notice is that it's broken up into blocks, with each block starting with two integers on a line, followed by a list of floating point numbers. The two integers are the block number, and the number of coefficients the block contains. The number of coefficients, is the same as NCOEFF from the header file. And remember that the number of coefficients isn't always equal to the actual count of numbers. For example the DE405 blocks are always padded with two zeros, so there are actually 1020 numbers in each block. Other ephemeris versions will have different padding standards.

The first two numbers of each block (for every ephemeris version) are the Julian Days for which the block is valid. For example, the first block in the ascp2020.405 file is valid from 2458832.5 to 2458864.5.

The not-so obvious breakdown of the coefficients starts here. This is where you need the data from "GROUP 1050" given in the header file, as well as the knowledge of how many components each series has (given in Table 1 above). Using Mercury as an example, it has a starting offset of "3" (line 1 of GROUP 1050), so its coefficients start with -0.468225142464447618D+08.

Line 2 of GROUP 1050 shows how many coefficients each component has, for Mercury that's 14. We know from Table 1 that Mercury has three components: X, Y and Z.

Using the first block of ascp2020.405 as an example, for the first subinterval (explained in the next section) the Mercury's coefficients for X, Y, and Z are:

X:-0.468225142464447618D+08 0.855287673857185431D+07 0.612484375662173959D+06 -0.271197032404459242D+05 0.175867546549079435D+03 -0.162623577504849326D+01 -0.763823621681322784D+00 0.441800821623631670D-01 -0.225469494234295390D-02 0.850807253766397409D-04 -0.218587309569990039D-05 0.850227664778876227D-09 0.490853028896027249D-08 -0.398938554977480726D-09

Y: -0.427358720430963337D+08 -0.935967597832447290D+07 0.576081445763668744D+06 0.126923156433521799D+05 -0.614832892580008775D+03 0.253491067956385763D+02 -0.982133401189059008D+00 0.199660889047113960D-01 -0.651480991977301560D-04 -0.458218756391133968D-04 0.341552453779017768D-05 -0.189986098666370616D-06 0.790647061064605362D-08 -0.238712559861458182D-09

Z: -0.181325881156450957D+08 -0.588670663881796692D+07 0.244250146726330160D+06 0.959132849803752833D+04 -0.346669478216073912D+03 0.137098683499955740D+02 -0.445471694283314124D+00 0.608610133598603501D-02 0.198916433185458984D-03 -0.332970361189953508D-04 0.205113183990416121D-05 -0.101577300506667819D-06 0.371477546608301747D-08 -0.861673297708620472D-10

Line 3 of GROUP 1050 is the number of subintervals a series has in each block. As explained in the "Valid Dates" section, each block has data for only a specific period of time (usually 32 days). Some of the planets, like Mercury, have to be divided into even smaller sections. We see that from line 3, in the first colum of GROUP 1050, that Mercury is divided into 4 subintervals. This means that all of the coefficients for all of the components are included 4 times, each set valid only for 1/4th the interval in which the entire block is valid for.

Using Mercury and the first block of ascp2020.405 as an example. It has 4 subintervals, the block is valid for 32 days, so each subinterval for Mercury is valid for 8 days. So, 2458832.5 is the start of the first subinterval, and 2458840.5 is the start of the second subinterval, and so on. So, if we were asked to compute the position of Mercury on JD=2458850.5, we would see that date falls in the third subinterval (subinterval 2 if counting the first as 0). We would then compute:

```
lengthOfSubinterval = daysPerBlock / numberOfSubintervals
lengthOfSubinterval = 32 / 4 = 8
subinterval = floor( (JD-blockStartDate)/lengthOfSubinterval )
subinterval = floor( (2458850.5 - 2458832.5)/8)
subinterval = floor( 2.25 )
subinterval = 2;
offset=subinterval*numberOfCoefficients*numberOfProperties+seriesStartOffset;
offset=2 * 14 * 3 + 3 = 87
```

So we would start with the 87th number in the list of floats, and take the X, Y, and Z coefficients from there. Remember that if you're using arrays, most languages start indexes with 0, so the 87th number is array index 86.

X: 0.230446411715880504D+04 0.133726736662702635D+08 -0.782187090879053358D+04 -0.267678745522568279D+05 -0.227070698075548364D+03 -0.142012340261296774D+02 -0.924872006275108544D-01 0.431659104815666252D-02 0.356917634561652571D-03 0.302564651657819373D-04 0.980701702776103911D-06 0.505819702568259545D-07 0.113034198242195379D-08 0.323800745882515925D-10

Y: -0.593914454531169161D+08 0.138391312173493067D+07 0.725419090211108793D+06 0.139471465250126903D+04 -0.290917422263861397D+03 -0.635064566332839320D+01 -0.646844700926034299D+00 -0.120797394835047579D-01 -0.681164244772722110D-03 -0.783160742259704191D-05 -0.953933699143903451D-07 0.170514411319974421D-07 0.132846579503924915D-08 0.625629348278007546D-10

Z: -0.318846685234071501D+08 -0.647159726192409638D+06 0.388325111500594590D+06 0.351975238047553557D+04 -0.131868705094903135D+03 -0.192040059987689204D+01 -0.335952534459033059D+00 -0.690035617434751804D-02 -0.400870301836372738D-03 -0.731991272299537233D-05 -0.152615685994738755D-06 0.386553675297770635D-08 0.592487943320094233D-09 0.300642066655273442D-10

All that work just to find a handful of numbers we want to multiply against. The good news is the JPL has already done all of the hard work from here, there's just a couple of more steps to get the position.

The Chebyshev polynomials don't use Julian Days for computation, instead they require a number from -1 to +1. -1 would represent the time of the beginning
of a subinterval, and +1 represents the end. So we have to scale our Julian Day based on the start and end times of the subinterval. It's important to
remember to use the times for which the **subinterval** is valid, not the entire block. So we must first compute the valid range for the subinterval.
Using the example of Mercury and JD=2458850.5 from above.

```
validStart= blockStart + subinterval * lengthOfSubinterval = 2458848.5
validEnd= validStart + lengthOfSubinterval = 2458856.5
```

In order to preserve maximum accuracy of the floating point JD, it is necessary to subtract off the validStart before normalization:

```
temp = JD - validStart
temp = 2458850.5 - 2458848.5 = 2.0
```

You may have noticed that many programs use a two part JD in order to increase the accuracy of the JD. In all previous instances where we used the JD, you can just add the two parts together without loss of accuracy. It is here at this step where it makes a difference.

```
[alternate two part JD method]
temp = JDpart1 - validStart + JDpart2
```

Now compute the normalized time variable which will be used in the Chebyshev polynomial:

```
x = temp / lengthOfSubinterval * 2.0 - 1.0
x = 2.0 / 8.0 * 2.0 - 1.0
x = -0.5
```

Computing the Chebyshev polynomial involves expanding the polynomial using the algorithm below. This is equation 14.20 from the Explanatory Supplement to the Astronomical Almanac, or alternatively presented on the Wikipedia Page:

```
T(0) = 1
T(1) = x
T(n) = 2 * x * T(n-1) - T(n-2)
```

Where n is the number of coefficients for the property we're computing.

Continuing our example of computing the position of Mercury for JD=2458850.5. We have x = -0.5 from above, and 14 coefficients, so n = 14. Note we're using array indexes starting at 0, so n will actually stop at 13.

T(0) = 1 T(1) = -0.5 T(2) = 2 * -0.5 * -0.5 - 1 T(3) = 2 * -0.5 * -0.5 - -0.5 T(4) = 2 * -0.5 * 1 - -0.5 T(5) = 2 * -0.5 * -0.5 - 1 T(6) = 2 * -0.5 * -0.5 - -0.5 T(7) = 2 * -0.5 * 1 - -0.5 T(8) = 2 * -0.5 * -0.5 - 1 T(9) = 2 * -0.5 * -0.5 - -0.5 T(10) = 2 * -0.5 * 1 - -0.5 T(11) = 2 * -0.5 * -0.5 - 1 T(12) = 2 * -0.5 * -0.5 - -0.5 T(13) = 2 * -0.5 * 1 - -0.5

Now it's a simple matter of multiplying each coefficient by it's t[n] counterpart and summing them together. We loop through the coefficients in reverse order (from smallest to largest) to avoid rounding errors. The end routine that does all of the work is simply (in JavaScript):

```
function computePolynomial(x,coefficients){
let T=new Array();
T[0]=1;
T[1]=x;
for(let n=2;n<14;n++) {
T[n]=2*x*T[n-1] - T[n-2];
}
let v=0;
for(let i=coefficients.length-1;i>=0;i--){
v+=T[i]*coefficients[i];
}
return v;
}
```

A function which calls that routine with our example data is:

```
function computeExamplePolynomials(){
let X=[0.230446411715880504E+04, 0.133726736662702635E+08, -0.782187090879053358E+04, -0.267678745522568279E+05,
-0.227070698075548364E+03, -0.142012340261296774E+02, -0.924872006275108544E-01, 0.431659104815666252E-02,
0.356917634561652571E-03, 0.302564651657819373E-04, 0.980701702776103911E-06, 0.505819702568259545E-07,
0.113034198242195379E-08, 0.323800745882515925E-10];
let Y=[-0.593914454531169161E+08, 0.138391312173493067E+07, 0.725419090211108793E+06, 0.139471465250126903E+04,
-0.290917422263861397E+03, -0.635064566332839320E+01, -0.646844700926034299E+00, -0.120797394835047579E-01,
-0.681164244772722110E-03, -0.783160742259704191E-05, -0.953933699143903451E-07, 0.170514411319974421E-07,
0.132846579503924915E-08, 0.625629348278007546E-10];
let Z=[-0.318846685234071501E+08, -0.647159726192409638E+06, 0.388325111500594590E+06, 0.351975238047553557E+04,
-0.131868705094903135E+03, -0.192040059987689204E+01, -0.335952534459033059E+00, -0.690035617434751804E-02,
-0.400870301836372738E-03, -0.731991272299537233E-05, -0.152615685994738755E-06, 0.386553675297770635E-08,
0.592487943320094233E-09, 0.300642066655273442E-10];
let x=-0.5;
console.log(computePolynomial(x,X));
console.log(computePolynomial(x,Y));
console.log(computePolynomial(x,Z));
}
```

And the final output for our example:

```
X = -6706768.766943997 km
Y = -60444568.85087551 km
Z = -31751664.901437085 km
```

You thought we were done? In many cases it's required to also compute the velocity, so that relativistic effects on the observer's time can be accounted for. The ephemeris data doesn't contain any data for velocity, but, fortunately velocity is just the first derivitive of position. So all we need to do is compute the derivitive. I won't list out all of the polynomials again, but the final function which computes both position and velocity, as well as a slightly modified calling function to finish the derivitive out:

```
function computePolynomial(x,coefficients){
//Equation 14.20 from Explanetory Supplement 3rd ed.
const t=new Array();
t[0]=1;
t[1]=x;
for(let n=2;n
```=0;i--){
position+=coefficients[i]*t[i];
}
//Compute velocity (just the derivitave of the above)
const v=new Array();
v[0]=0;
v[1]=1;
v[2]=4*x;
for(let n=3;n=0;i--){
velocity+=v[i]*coefficients[i];
}
let retval=new Array();
retval[0]=position;
retval[1]=velocity;
return retval;
}
function computeExamplePolynomials2(){
let X=[0.230446411715880504E+04, 0.133726736662702635E+08, -0.782187090879053358E+04, -0.267678745522568279E+05,
-0.227070698075548364E+03, -0.142012340261296774E+02, -0.924872006275108544E-01, 0.431659104815666252E-02,
0.356917634561652571E-03, 0.302564651657819373E-04, 0.980701702776103911E-06, 0.505819702568259545E-07,
0.113034198242195379E-08, 0.323800745882515925E-10];
let Y=[-0.593914454531169161E+08, 0.138391312173493067E+07, 0.725419090211108793E+06, 0.139471465250126903E+04,
-0.290917422263861397E+03, -0.635064566332839320E+01, -0.646844700926034299E+00, -0.120797394835047579E-01,
-0.681164244772722110E-03, -0.783160742259704191E-05, -0.953933699143903451E-07, 0.170514411319974421E-07,
0.132846579503924915E-08, 0.625629348278007546E-10];
let Z=[-0.318846685234071501E+08, -0.647159726192409638E+06, 0.388325111500594590E+06, 0.351975238047553557E+04,
-0.131868705094903135E+03, -0.192040059987689204E+01, -0.335952534459033059E+00, -0.690035617434751804E-02,
-0.400870301836372738E-03, -0.731991272299537233E-05, -0.152615685994738755E-06, 0.386553675297770635E-08,
0.592487943320094233E-09, 0.300642066655273442E-10];
let x=-0.5;
let temp;
temp=computePolynomial2(x,X);
temp[1]=temp[1]*(2.0*4.0/32.0);
console.log(temp);
temp=computePolynomial2(x,Y);
temp[1]=temp[1]*(2.0*4.0/32.0);
console.log(temp);
temp=computePolynomial2(x,Z);
temp[1]=temp[1]*(2.0*4.0/32.0);
console.log(temp);
}

And the final output with velocities is:

X = -6706768.766943997 km dX = 3346870.03970893 km/day Y = -60444568.85087551 km dY = -17014.263564507186 km/day Z = -31751664.901437085 km dZ = -356081.96677701955 km/day

And these match to within 1km of the values from JPL Horizons.

Most ephemeris implementations output data in AU and AU/Day, while the above works in kilometers and km/day which is what the JPL Ephemeris uses internally. If you want to match the output of other programs, you will need to covert.

The code produced here was written primarily for clarity. There are many optimizations that can be done to make it faster, though you may find that it works well enough as is. If you're looking for places for optimization, the first place to look is the generation of the arrays in "computePolynomial". The arrays are regenerated for every call, but it is highly likely that you'll make multiple calls to this function with the same value of "x", so it would help quite a bit to cache the array and check to see if it needs to be recomputed or not. And, if you always generate the maxiumim entries required by any series, you can reuse the array for every planet/series.

Another good place to look for optimization is "computePolynomial" accepts an array of polynomials. Which requires copying of data from the main data array into a temporary array. It would be much faster to just pass an index into the main data array.

The JPL ASCII data files are available on the JPL FTP Site

Not all data is available for every DE version. Some versions specify a larger interval than there are files on the JPL FTP site. But there are test cases for these unavailable data ranges in the testpo.xxx files. So you will encounter errors if you blindly run all tests.