12
$\begingroup$

I'm attempting to do some basic planet position calculations using the JPL ephemeris data. The first step is to properly parse the files. Here's my understanding so far:

The header file gives a matrix, for example:

 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 documentation says that this can be interpreted as follows:

Word (1,i) is the starting location in each data record of the chebychev coefficients belonging to the ith item. Word (2,i) is the number of chebychev coefficients per component of the ith item, and Word (3,i) is the number of complete sets of coefficients in each data record for the ith item.

And later on:

The first two double precision words in each data record contain the Julian date of earliest data in record and the Julian date of latest data in record.

The remaining data are chebychev position coefficients for each component of each body on the tape. The chebychev coefficients for the planets represent the solar system barycentric positions of the centers of the planetary systems.

There are three Cartesian components (x, y, z), for each of the items #1-11

So each planet has $N$ coefficients in each of $n$ subintervals for each of the 3 coordinates ($N = 14$ and $n=4$ for the first planet in the data above).

My question now is, what is the order of the coefficients in each subinterval? For example if $X_i^k, Y_i^k$ and $Z_i^k$ are $k$th coefficient for $(x,y,z)$ in the $i$th subinterval, the coefficients could be logically ordered in any of the following ways (and probably others):

  • $X_1^1, ..., X_1^N, Y_1^1, Y_1^2, ..., Y_1^N, Z_1^1, ..., Z_1^N, X_2^1, ...$
  • $X_1^1, ..., X_1^N, X_2^1, ... X_n^N, Y^1_1, ..., Y_n^N, Z^1_1, ..., Z_n^N$
  • $X_1^1, Y_1^1, Z_1^1, X_1^2, Y_1^2, Z_1^2, ...$

Is there any documentation that goes over this order?

$\endgroup$
0

4 Answers 4

14
$\begingroup$

You can parse and use those ascii JPL ephemeris files.

Warning: This is an exercise for people who like to torture themselves. If you aren't into self-flagellation, you should do what Mark Adler wrote and use the SPICE toolkit. I went through this self-flagellation over a decade ago, when SPICE was closed and wasn't what it is now. Now people are thinking of killing a thousand lines of code or so by replacing my work with SPICE. And that is exactly what they should do.


Consider the file ascp2200.405. (I chose that file because it's relatively small.) The first line is

     1  1018

The "1" denotes the block number within the file. The "1018" denotes the number of numbers in the block. The next 340 lines contain those 1018 numbers (plus two more for padding). The first two numbers represent the start and end Julian dates of the block, in JPL's own ephemeris time. If you don't care about millisecond-ish errors, you can use TAI in lieu of JPL's $T_{\text{eph}}$.

The first 168 numbers in that block pertain to Mercury, the next 60 to Venus, the next 78 to the Earth-Moon barycenter, and so on. This is specified in the GROUP 1050 section of the header file. You can deduce the dimensionality of each item from that section. Alternatively, you should know that position is three-dimensional, while librations are two dimensional.

The 168 numbers that describe the position of Mercury are organized in four sub-blocks, each comprising 3*14=42 numbers. (Regarding those magic numbers: 14 is from the GROUP 1050 data for Mercury, and 3 is the dimensionality of the position data.) Of those 42 numbers, the first 14 are the coefficients for the 13th order Chebyshev polynomial coefficients that describes the x position, the next 14 describe the y position, and the last 14 describe the z position.

Those 14 numbers are Chebyshev polynomial coefficients with time (properly scaled) as the independent variable. Chebyshev polynomials are functions limited to the domain [-1,1]. To find the position of Mercury at some time $t$, you'll first need to find the block that pertains to that time. Then you'll need to find the sub-block (Mercury has four sets of coefficients) that pertains to that time. Next you'll have to offset the time from the middle of that sub-block and scale the result so that -1 represents the start of the sub-block and 1 represents the end. Finally, you'll need to compute $\sum_{k=0}^{13} C_n T_n(t)$ for each of the three elements of the position vector, where $C_n$ are the Chebyshev polynomial coefficients and the $T_n(t)$ are the Chebyshev polynomials: $T_0(t)=1$, $T_1(t)=t$, and $T_{n+1}(t) = 2tT_n(t)-T_{n-1}(t)$. The results are in kilometers, except for items pertaining to rotation. This may or may not be what you want.

To make matters even more fun (self-flagellation fun, that is), what if you want to know the position of Mars relative to the Earth? The JPL ephemeris file doesn't have coefficients for the Earth. It has coefficients for the Earth-Moon barycenter (relative to the solar system barycenter) and for the Moon (relative to the center of the Earth). You'll need to calculate the position of the Earth by calculating the position of Mars, the position of the Earth-Moon barycenter, and the position of the Moon relative to the Earth. Next you'll need to calculate the position of the Earth (relative to the solar system barycenter) using the Earth-Moon mass ratio (a magic number that's in the header file). With this, you can finally compute the position of Mars relative to the Earth. That's the instantaneous relative position. If you care about where you would see Mars in the sky, you'll need to worry about where Mars was at the time the light reaches you on the Earth.


Alternatively, you could follow Mark Adler's advice: Use the SPICE toolkit.

$\endgroup$
0
13
$\begingroup$

No, the first step is to download the SPICE Toolkit. It will read the ephemeris and generate the data you need for you, as well as do much, much more. There are many formats used in SPICE kernel files, not just the particular one you're looking at, which is probably a Chebyshev position-only array (type 2). There are 16 types. The values can be stored little or big endian. There can be large comment sections in the file. There are several bodies in each file. The bodies will have different degrees of resolution in time. It will take some work for you to decode the thing yourself, which is just wasted time since the toolkit will do all that for you.

Download the toolkit.

To answer to your question, if you have successfully found a single type 2 Chebyshev position-only array record and have arrived at the polynomial coefficients (having read MID and RADIUS for that record), is that: the sequence of doubles starts with the constant term of the polynomial for $x$, and continues with the linear coefficient of $t$, until that polynomial is complete, followed by the same thing for $y$ and $z$. The doubles are IEEE 64-bit floating point numbers in little or big-endian order depending on what it says in the header.

$t$ is ephemeris time in seconds, and the results are in km.

The format is documented here.

However you should download the toolkit.

$\endgroup$
5
  • $\begingroup$ The toolkit is a good idea, however I should have mentioned in my question that this is for my own personal interest, so I'd like to code up a rough version on my own. The format in the link you provided looks slightly different than the format described here. For one thing the time are provided in start and end, not mid and radius. Also the coefficients should be coefficients for Chebyshev polynomials, not the polynomial you described. $\endgroup$ Oct 30, 2015 at 16:31
  • $\begingroup$ The first two Chebyshev polynomials happen to be $1$ and $x$, so what I said is consistent with that. $\endgroup$
    – Mark Adler
    Oct 30, 2015 at 16:55
  • $\begingroup$ You seem to be conflating two different parts of the format. The meat of the data is a sequence of records where each are each mid, radius, x coeffs, y, coeffs, z coeffs. That sequence of records is followed by four doubles that are the initial epoch, interval length, elements in each record (from you which can compute the degree of the polynomials), and the number of records. $\endgroup$
    – Mark Adler
    Oct 30, 2015 at 16:58
  • $\begingroup$ Or maybe you are looking at some different sort of files. You should be using the SPK Ephemeris files such as de430.bsp or de431.bsp. $\endgroup$
    – Mark Adler
    Oct 30, 2015 at 16:59
  • 1
    $\begingroup$ The deep innards of the JPL ephemerides are Chebyshev polynomials of the first kind, not second: $T_0(x) = 1$, $T_1(x)=x$, $T_{n+1}(x) = 2xT_n(x)-T_{n-1}(x)$. $\endgroup$ Oct 30, 2015 at 18:13
3
$\begingroup$

The exact format of the SPICE Toolkit SPK files (also called BSP files) is called DAF (Double-precision Array File). It is a binary format made up of consecutive blocks of 1,024 bytes. The first record is always the file record, then it has some blocks serving as a comment area, then it has the array records, and finally a collection of data records.

[FILE RECORD] (1024 bytes) [FIRST COMMENT BLOCK] (1024 bytes) ... [LAST COMMENT BLOCK] (1024 bytes) [FIRST BLOCK OF SUMMARY RECORDS] (1024 bytes) ... [LAST BLOCK OF SUMMARY RECORDS] (1024 bytes) [FIRST BLOCK OF NAME RECORDS] (1024 bytes) ... [LAST BLOCK OF NAME RECORDS] (1024 bytes) [FIRST BLOCK OF ELEMENT RECORDS] (1024 bytes) ... [LAST BLOCK OF ELEMENT RECORDS] (1024 bytes)

All binary SPICE data is stored in DAF format (not only ephemerides). The specific kind of data (that is, the semantics of the records) depend on the info stored in the summary records and the element records. You need to know what kind of data you are reading to correctly interpret the records.

There are many ephemerides types. Planets are often Type II; satellites are often Type III; spacecraft are often Type I. Types II and III represent Chebyshev polynomials, and their difference is whether they use the derivative of the Chebyshev polynomial to calculate velocity (Type II), or whether a separate polynomial is used to calculate velocity (Type III).

In other words: a Type II SPK file will provide you with one polynomial of certain order whose value represents position, and whose derivative represents speed. A Type III SPK file will provide you with one polynomial for position and a different polynomial for velocity. Both types are Chebyshev polynomials.

For SPK files, the summary records will tell you where every record begins and ends, the ephemerides type, the timespan, the frame, the target, and the observer. Each specific records (which will most likely contain more than one block) look as follows:

+---------------+ | Record 1 | +---------------+ | Record 2 | +---------------+ . . . +---------------+ | Record N | +---------------+ | INIT | +---------------+ | INTLEN | +---------------+ | RSIZE | +---------------+ | N | +---------------+

and each record looks like

+------------------+ | MID | +------------------+ | RADIUS | +------------------+ | X coefficients | +------------------+ | Y coefficients | +------------------+ | Z coefficients | +------------------+

for Type II and like

+------------------+ | MID | +------------------+ | RADIUS | +------------------+ | X coefficients | +------------------+ | Y coefficients | +------------------+ | Z coefficients | +------------------+ | X' coefficients | +------------------+ | Y' coefficients | +------------------+ | Z' coefficients | +------------------+

for Type III.

The degree of the polynomial is ( RSIZE - 2 ) / 3 - 1. The MID tells you the interval mid-point (ephemerides time), and RADIUS tells you half the size of the interval. So its time span is [MID - RADIUS, MID + RADIUS].

For obtaining position you simply evaluate the Chebyshev at the corresponding time. For velocity you either evaluate the derivative or a different polynomial. Remember to normalize time before evaluating and de-normalize time after evaluating (this has gotten me many times).

Dealing directly with DAF can be a rewarding experience (for example: I have a thread-safe high-performance ephemerides evaluator), but it is easy to get wrong. And after a lot of work you only have a position or velocity in a given frame that you still need to rotate/translate/etc (for which you either need to re-invent SPICE or rely on it).

Unless you have a very good reason not to, just stick to SPICE.

A simple Python script

I just made a Python script that could help someone get started on the DAF binary format. It will simply dump the file record information and the contents (array summaries) contained in a DAF.

import mmap
import struct

RECLEN = 1024

def peek_spk(mem):
    # file record is always the first one

    # string data
    locidw = mem[0:8]
    locifn = mem[16:16 + 60]
    locfmt = mem[88:88+8]
    ftpstr = mem[699:699+28]

    # endianness
    fmt     = "<" if locfmt == "LTL-IEEE" else ">"
    int_fmt = fmt + "I"

    nd,    = struct.unpack_from(int_fmt, mem, offset=8)
    ni,    = struct.unpack_from(int_fmt, mem, offset=12)
    fward, = struct.unpack_from(int_fmt, mem, offset=76)
    bward, = struct.unpack_from(int_fmt, mem, offset=80)
    free,  = struct.unpack_from(int_fmt, mem, offset=84)

    print "locidw {0}".format(locidw)
    print "nd     {0}".format(nd)
    print "ni     {0}".format(ni)
    print "locifn {0}".format(locifn)
    print "fward  {0}".format(fward)
    print "bward  {0}".format(bward)
    print "free   {0}".format(free)
    print "locfmt {0}".format(locfmt)
    print "ftpstr {0}".format(repr(ftpstr))

    # let's do the first summary record only... we need to read nd
    # doubles and ni integers starting at offset fward

    offset = (fward - 1) * RECLEN
    sum_fmt = fmt + nd * "d" + ni * "I"
    size = nd + (ni + 1) / 2 # integer division
    while True:
        nxt, prv, nsum = struct.unpack_from(fmt + "ddd", mem, offset=offset)
        offset += 24 # skip three doubles
        for n in range(int(nsum)):
            print struct.unpack_from(sum_fmt, mem, offset = offset)
            offset += size * 8
        if nxt == 0:
            break


def peek(path):

    print "peeking into {0}".format(path)

    with open(path, "rb") as fp:
        mem = mmap.mmap(fp.fileno(), 0, access=mmap.PROT_READ)

    peek_spk(mem)


if __name__ == "__main__":
    import sys
    for path in sys.argv[1:]:
        peek(path)

For example: running it for jup310.bsp gives

peeking into /data/spice/jup310.bsp locidw DAF/SPK nd 2 ni 6 locifn NIO2SPK
fward 6 bward 6 free 122077514 locfmt LTL-IEEE ftpstr 'FTPSTR:\r:\n:\r\n:\r\x00:\x81:\x10\xce:ENDFTP' (-3155716758.8160305, 3155716867.183885, 501, 5, 1, 3, 897, 7208500) (-3155716758.8160305, 3155716867.183885, 502, 5, 1, 3, 7208501, 14367404) (-3155716758.8160305, 3155716867.183885, 503, 5, 1, 3, 14367405, 19140106) (-3155716758.8160305, 3155716867.183885, 504, 5, 1, 3, 19140107, 22451778) (-3155716758.8160305, 3155716867.183885, 505, 5, 1, 3, 22451779, 44074360) (-3155716758.8160305, 3155716867.183885, 514, 5, 1, 3, 44074361, 65696942) (-3155716758.8160305, 3155716867.183885, 515, 5, 1, 3, 65696943, 90825888) (-3155716758.8160305, 3155716867.183885, 516, 5, 1, 3, 90825889, 115954834) (-3155716758.8160305, 3155716867.183885, 599, 5, 1, 3, 115954835, 120922238) (-3155716758.8160305, 3155716867.183885, 3, 0, 1, 2, 120922239, 121109489) (-3155716758.8160305, 3155716867.183885, 5, 0, 1, 2, 121109490, 121168877) (-3155716758.8160305, 3155716867.183885, 10, 0, 1, 2, 121168878, 121328726) (-3155716758.8160305, 3155716867.183885, 399, 3, 1, 2, 121328727, 122077513)

That means that jup310 contains records spanning ET (-3155716758.8160305, 3155716867.183885) for some moons of Jupiter (501 through 519) relative to Jupiter's Barycenter (5), in J2000 format (1), and are type 3. Other records are for the Barycenters of Earth, Jupiter, and the Sun (3, 5, and 10) relative to the Solar System Barycenter (0) in J200 (1), and are type 2.

Each record tells you where in the file does each array begin and end. You can extend the script to skip to each location and read the entire record for later interpolation.

Available Software

I quickly created a simple implementation in Python which can read and work directly with Types II and III. It compares well to CSPICE.

https://gist.github.com/arrieta/c2b56f1e2277a6fede6d1afbc85095fb

$\endgroup$
1
  • $\begingroup$ this is a comprehensive and extremely useful answer. Thank you. I wonder if you could elaborate on the use of the factor 2 during polynomial evaluation i.e. why the input to the Horner method is 2x rather than x? $\endgroup$ May 3, 2018 at 17:59
1
$\begingroup$

I'm going to disagree with all of the previous answers. I would argue that implementing your own Ephemeris program is borderline trivial once you know the format of the files and have the relatively simple equations for using the polynomials for computing the positions. The task is well within the realm of a beginner programmer. The format of the file is the hardest part, only because there are a lot of pieces to put together, but not exactly difficult to do.

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 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
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.