What is the mathematics needed for understanding DE440 and translating it into xyz coordinates?

Question

I would like to write my own code to process NASA's DE440/DE441 ephemerides, but I do not know the mathematics needed to do this.

I understand packages such as 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?

Answer 1

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.