<p>This is largely a copy of my answer from <a href="https://space.stackexchange.com/questions/12506/what-is-the-exact-format-of-the-jpl-ephemeris-files">What is the exact format of the JPL ephemeris files?</a> with a few modifications.  The below is a description of how to use the ASCII files, which is the easiest to use.  The article linked below explains the binary format too, but knowledge of the ASCII format is required.  Once you have the correct coefficients, using them to compute the coordinates is quite trivial.</p>
<p>The article <a href="http://astrogreg.com/jpl-ephemeris-format/jpl-ephemeris-format.html" rel="nofollow noreferrer">Format of the JPL Ephemeris Files</a> has a pretty thorough walk through, with examples, and a <a href="http://astrogreg.com/jpl-ephemeris-format/de405javascript.html" rel="nofollow noreferrer">JavaScript implementation</a> as an example.  Note the whole class that does the coefficient lookups and computations is less than 100 lines of code.  Additionally there are example implementations in the Github repository <a href="https://github.com/gmiller123456/jpl-development-ephemeris" rel="nofollow noreferrer">gmiller123456/jpl-development-ephemeris</a> which currently has implementations in Python, Perl, C# and JavaScript.</p>
<p>Like you've already figured out, the &quot;GROUP 1050&quot; section of the header contains a lot of the keys to the format of the file.  This is the section for DE405</p>
<pre><code>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
</code></pre>
<p>The order, and number of components which these correspond to is given in the <a href="ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/ascii_format.txt" rel="nofollow noreferrer">Ascii Format</a> document from JPL.  Each column represents the values for a given planet (or other property):</p>
<pre><code>#    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)
</code></pre>
<p>Each of the <a href="ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/" rel="nofollow noreferrer">ASCII Files</a> is divided up into blocks, valid for the number of days given in the header file.  The first two numbers in the block are the Julian dates for which that block is valid, so finding the corresponding block to a given Julian date is trivial.</p>
<p>From GROUP 1050 in the header, the numbers in each column correspond to:</p>
<pre><code>1. The start offset in the block (starting at 1)
2. Number of coefficients for each property
3. The number of subintervals
</code></pre>
<p>So, to compute the offset into the block for a given planet:</p>
<pre><code>lengthOfSubinterval = daysPerBlock / numberOfSubintervals
subinterval = floor( (JD-blockStartDate)/lengthOfSubinterval )
offset=subinterval*numberOfCoefficients*numberOfProperties+seriesStartOffset;
</code></pre>
<p>So, for planets, you'll have all of the X coefficients, then all of the Y coefficients, and finally all of the Z coefficients.  For example the coefficients for Mercury for JD=2458850.5 are given in the code below.</p>
<p>To expand these coefficients into expand them into polynomials using equations 14.20 from the Explanatory Supplement to the Astronomical Almanac:
<span class="math-container">$$
T_0(x)=1 \\
T_1(x)=x \\
T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)
$$</span>
where <span class="math-container">$x$</span> is the time normalized from -1.0 to 1.0:</p>
<pre><code>validStart= blockStart + subinterval * lengthOfSubinterval
temp = JD - validStart
x = temp / lengthOfSubinterval * 2.0 - 1.0
</code></pre>
<p>And then multiply each T_n(x) against its respective coefficient.  It may sound complicated mathematically, but a code example should make it quite clear.  To borrow the example from the <a href="http://astrogreg.com/jpl-ephemeris-format/jpl-ephemeris-format.html" rel="nofollow noreferrer">Format of the JPL Ephemeris Files</a> article, this actually computes the positions:</p>
<pre><code>function computePolynomial(x,coefficients){
  let T=new Array();

  T[0]=1;
  T[1]=x;
  for(let n=2;n&lt;14;n++)  {
    T[n]=2*x*T[n-1] - T[n-2];
  }

  let v=0;
  for(let i=coefficients.length-1;i&gt;=0;i--){
    v+=T[i]*coefficients[i];
  }
  return v;
}    

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));
}
</code></pre>
<p>Calling the computeExamplePolynomials() function above, produces the values below:</p>
<pre><code>X = -6706768.766943997 km
Y = -60444568.85087551 km
Z = -31751664.901437085 km
</code></pre>
<p>The article also includes information on to compute the velocity of each planet, which is required for stellar aberration corrections.</p>

CC BY-SA 4.0
2023-06-27T12:16:47.070Z
2023-06-27T12:36:25.403Z
2023-06-27T12:36:25.403Z
Greg Miller
<p>I would like to write my own code to process <a href="https://ssd.jpl.nasa.gov/doc/de440_de441.html" rel="nofollow noreferrer">NASA's DE440/DE441 ephemerides</a>, but I do not know the mathematics needed to do this.</p>
<p>I understand packages such as <a href="https://rhodesmill.org/pyephem/" rel="nofollow noreferrer">pyephem</a> have implemented variations of this, but I'd like to understand how to do this myself from the mathematics involved.</p>
<p>Is there a good reference, or can someone please provide a guide on how I can proceed?</p>

CC BY-SA 4.0
2023-06-27T05:13:42.870Z
Some excellent reading can be found by [searching for "Chebychev" in Space Exploration SE](https://space.stackexchange.com/search?q=chebychev). The DE's and JPL's SPICE kernels use a specific custom kind of Chebychev polynomial that minimizes interpolation error between time points in the [epeherides](https://english.stackexchange.com/q/370655/217285).
uhoh
CC BY-SA 4.0
2023-06-27T05:15:17.767Z
You could for example start with the open source Python package [Skyfield](https://rhodesmill.org/skyfield/) which uses and manages the JPL DE's directly and work backwards, or you might find Skyfield so darn easy to use you think twice about reinventing the wheel. There's a lot of work that has to go into understanding how to manage the way the coefficients are stored in discreet chunks, a lot of non-mathematical, uninteresting housekeeping to do it exactly correctly. BUT has a fun project just to do it, I think it will be quite fun and rewarding. *Kudos to you!*
uhoh
CC BY-SA 4.0
2023-06-27T05:19:55.517Z
I added the [tag:skyfield] and [tag:python] tags because we don't have a Pyephem tag. Both Pyephem and Skyfield are maintained by the same person, who recommends that if at all possible, use Skyfield - Pyephem is maintained but no longer being developed as far as I understand, and it does not use the DE's. It uses something older. [PyEphem under the hood - how does it calculate position of planets?](https://stackoverflow.com/q/34487170)
uhoh
CC BY-SA 4.0
2023-06-27T05:23:13.653Z
Just for example, see answers to [Is there a way to extract the Chebyshev coefficients for a body from a SPICE kernel?](https://space.stackexchange.com/q/10490/12102) and also I've just asked: [What is the difference between a "SPICE kernel" and a JPL Development Ephemeris? Is a DE just a bunch of kernels on a string?](https://astronomy.stackexchange.com/q/54110/7982)
uhoh
CC BY-SA 4.0
2023-06-27T11:42:39.337Z
This article explains the format and the algorithm needed to produce the coordinates, and has sample code in a few languages: https://www.celestialprogramming.com/jpl-ephemeris-format/jpl-ephemeris-format.html
Greg Miller
<orbital-mechanics><mathematics><python><celestial-mechanics><skyfield>
What is the mathematics needed for understanding DE440 and translating it into xyz coordinates?
Ana
# What is the mathematics needed for understanding DE440 and translating it into xyz coordinates?

(Asked by: username_0 on 2023-06-27)

I would like to write my own code to process [NASA's DE440/DE441 ephemerides](https://ssd.jpl.nasa.gov/doc/de440_de441.html), but I do not know the mathematics needed to do this.

I understand packages such as [pyephem](https://rhodesmill.org/pyephem/) have implemented variations of this, but I'd like to understand how to do this myself from the mathematics involved.

Is there a good reference, or can someone please provide a guide on how I can proceed?

**Tags:** orbital-mechanics, mathematics, python, celestial-mechanics, skyfield
> **Comments:**
> **username_1**: Some excellent reading can be found by [searching for "Chebychev" in Space Exploration SE](https://space.stackexchange.com/search?q=chebychev). The DE's and JPL's SPICE kernels use a specific custom kind of Chebychev polynomial that minimizes interpolation error between time points in the [epeherides](https://english.stackexchange.com/q/370655/217285).  
> **username_1**: You could for example start with the open source Python package [Skyfield](https://rhodesmill.org/skyfield/) which uses and manages the JPL DE's directly and work backwards, or you might find Skyfield so darn easy to use you think twice about reinventing the wheel. There's a lot of work that has to go into understanding how to manage the way the coefficients are stored in discreet chunks, a lot of non-mathematical, uninteresting housekeeping to do it exactly correctly. BUT has a fun project just to do it, I think it will be quite fun and rewarding. *Kudos to you!*  
> **username_1**: I added the [tag:skyfield] and [tag:python] tags because we don't have a Pyephem tag. Both Pyephem and Skyfield are maintained by the same person, who recommends that if at all possible, use Skyfield - Pyephem is maintained but no longer being developed as far as I understand, and it does not use the DE's. It uses something older. [PyEphem under the hood - how does it calculate position of planets?](https://stackoverflow.com/q/34487170)  
> **username_1**: Just for example, see answers to [Is there a way to extract the Chebyshev coefficients for a body from a SPICE kernel?](https://space.stackexchange.com/q/10490/12102) and also I've just asked: [What is the difference between a "SPICE kernel" and a JPL Development Ephemeris? Is a DE just a bunch of kernels on a string?](https://astronomy.stackexchange.com/q/54110/7982)  
> **username_2**: This article explains the format and the algorithm needed to produce the coordinates, and has sample code in a few languages: https://www.celestialprogramming.com/jpl-ephemeris-format/jpl-ephemeris-format.html  

---

## 1 answers

### Answer by username_2 on 2023-06-27. Score: 3

This is largely a copy of my answer from What is the exact format of the JPL ephemeris files? with a few modifications.  The below is a description of how to use the ASCII files, which is the easiest to use.  The article linked below explains the binary format too, but knowledge of the ASCII format is required.  Once you have the correct coefficients, using them to compute the coordinates is quite trivial.
The article Format of the JPL Ephemeris Files has a pretty thorough walk through, with examples, and a JavaScript implementation as an example.  Note the whole class that does the coefficient lookups and computations is less than 100 lines of code.  Additionally there are example implementations in the Github repository gmiller123456/jpl-development-ephemeris which currently has implementations in Python, Perl, C# and JavaScript.
Like you've already figured out, the "GROUP 1050" section of the header contains a lot of the keys to the format of the file.  This is the section for DE405
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

The order, and number of components which these correspond to is given in the Ascii Format document from JPL.  Each column represents the values for a given planet (or other property):
#    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)

Each of the ASCII Files is divided up into blocks, valid for the number of days given in the header file.  The first two numbers in the block are the Julian dates for which that block is valid, so finding the corresponding block to a given Julian date is trivial.
From GROUP 1050 in the header, the numbers in each column correspond to:
1. The start offset in the block (starting at 1)
2. Number of coefficients for each property
3. The number of subintervals

So, to compute the offset into the block for a given planet:
lengthOfSubinterval = daysPerBlock / numberOfSubintervals
subinterval = floor( (JD-blockStartDate)/lengthOfSubinterval )
offset=subinterval*numberOfCoefficients*numberOfProperties+seriesStartOffset;

So, for planets, you'll have all of the X coefficients, then all of the Y coefficients, and finally all of the Z coefficients.  For example the coefficients for Mercury for JD=2458850.5 are given in the code below.
To expand these coefficients into expand them into polynomials using equations 14.20 from the Explanatory Supplement to the Astronomical Almanac:
$$
T_0(x)=1 \\
T_1(x)=x \\
T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)
$$
where $x$ is the time normalized from -1.0 to 1.0:
validStart= blockStart + subinterval * lengthOfSubinterval
temp = JD - validStart
x = temp / lengthOfSubinterval * 2.0 - 1.0

And then multiply each T_n(x) against its respective coefficient.  It may sound complicated mathematically, but a code example should make it quite clear.  To borrow the example from the Format of the JPL Ephemeris Files article, this actually computes the positions:
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;
}    

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));
}

Calling the computeExamplePolynomials() function above, produces the values below:
X = -6706768.766943997 km
Y = -60444568.85087551 km
Z = -31751664.901437085 km

The article also includes information on to compute the velocity of each planet, which is required for stellar aberration corrections.


---

