Example demonstrating how could one work with SPICE SPK files in their original DAF format.

spk.py

"""spk.py

This is an example meant to demonstrate how could one work directly
with SPICE SPK files in their original DAF format.

It can correctly compute the position and velocity provided by Type II
and Type III ephemerides. As a basis for comparison it calculates the
states of the Galilean satellites for a relatively large number of
epochs and compares it agains CSPICE (which you need to have available
as a shared library if you want to run the test case, but you don't
need CSPICE at all to use the SPK class; it is standalone).

It is not meant to be high-performance. It is not meant for production
use. It is not mean to be well documented. It is not meant to be
cleanly implemented. It is just a quick example.

(C) 2017 Nabla Zero Labs <[email protected]>

Released under the terms of The MIT License:

Copyright 2017 Nabla Zero Labs

Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

"""

# Python Standard Library
import array
import math
import mmap
import struct

def chebyshev(order, x, data):
    """Evaluate a Chebyshev polynomial"""
    two_x = 2 * x
    bkp2 = data[order]
    bkp1 = two_x * bkp2 + data[order - 1]
    for n in xrange(order - 2, 0, -1):
        bk = data[n] + two_x * bkp1 - bkp2
        bkp2 = bkp1
        bkp1 = bk
    return data[0] + x * bkp1 - bkp2

def der_chebyshev(order, x, data):
    """Evaluate the derivative of a Chebyshev polynomial"""
    two_x = 2 * x
    bkp2 = order * data[order]
    bkp1 = two_x * bkp2 + (order - 1) * data[order - 1]
    for n in xrange(order - 2, 1, -1):
        bk = n * data[n] + two_x * bkp1 - bkp2
        bkp2 = bkp1
        bkp1 = bk
    return data[1] + two_x * bkp1 - bkp2

class SPK(object):

    RECLEN = 1024


    def __init__(self, path):
        self._path = path
        with open(path, "rb") as fp:
            self._mem =  mmap.mmap(fp.fileno(), 0, access=mmap.PROT_READ)

        self._prepare()

    def _prepare(self):

        if self._mem[0:7] != "DAF/SPK":
            msg = "{0} does not appear to be an SPK file"
            raise RuntimeError(msg.format(self._path))

        # deal with data endianness
        locfmt = self._mem[88:88+8]
        self._fmt = "<" if locfmt == "LTL-IEEE" else ">"

        int_fmt = self._fmt + "I"
        self._nd,    = struct.unpack_from(int_fmt, self._mem, offset=8)
        self._ni,    = struct.unpack_from(int_fmt, self._mem, offset=12)
        self._fward, = struct.unpack_from(int_fmt, self._mem, offset=76)
        self._bward, = struct.unpack_from(int_fmt, self._mem, offset=80)

        self._extract_summary_records()

    def _extract_summary_records(self):
        summary_offset = (self._fward - 1) * SPK.RECLEN
        summary_fmt = self._fmt + self._nd * "d" + self._ni * "I"
        summary_size = self._nd + (self._ni + 1) / 2 # integer division

        self._summary_records = []

        while True:
            nxt, prv, nsum = struct.unpack_from(self._fmt + "ddd", self._mem,
                                                offset=summary_offset)
            summary_offset += 24 # skip three doubles
            for n in range(int(nsum)):
                record = struct.unpack_from(summary_fmt, self._mem, offset=summary_offset)
                self._summary_records.append(record)
                summary_offset += summary_size * 8 # bytes
            if nxt == 0:
                break

        # TODO: add more information to the record (such as INIT,
        # INTLEN, RSIZE, and N) to avoid re-reading it every time.

    def state(self, et, target, observer):
        type, order, data = self.record(et, target, observer)
        if type == 2:
            return self.state_type_2(et, order, data)
        elif type == 3:
            return self.state_type_3(et, order, data)
        else:
            raise RuntimeError("Only Types II and III are implemented")

    def state_type_2(self, et, order, data):
        # todo: improve slicing... this is too stupid (but it works)
        tau = (et - data[0]) / data[1]
        beg = 2
        end = beg + order + 1
        deg = order  + 1
        factor = 1.0 / data[1] # unscale time dimension
        return (
            chebyshev(order, tau, data[2           : 2 + deg]),
            chebyshev(order, tau, data[2 + deg     : 2 + 2 * deg]),
            chebyshev(order, tau, data[2 + 2 * deg : 2 + 3 * deg]),
            # type 2 uses derivative on the same polynomial
            der_chebyshev(order, tau, data[2           : 2 + deg]) * factor,
            der_chebyshev(order, tau, data[2 + deg     : 2 + 2 * deg]) * factor,
            der_chebyshev(order, tau,  data[2 + 2 * deg : 2 + 3 * deg]) * factor,
            )

    def state_type_3(self, et, order, data):
        # todo: improve slicing... this is too stupid (but it works)
        tau = (et - data[0]) / data[1]
        beg = 2
        end = beg + order + 1
        deg = order  + 1
        return (
            chebyshev(order, tau, data[2           : 2 + deg]),
            chebyshev(order, tau, data[2 + deg     : 2 + 2 * deg]),
            chebyshev(order, tau, data[2 + 2 * deg : 2 + 3 * deg]),
            chebyshev(order, tau, data[2 + 3 * deg : 2 + 4 * deg]),
            chebyshev(order, tau, data[2 + 4 * deg : 2 + 5 * deg]),
            chebyshev(order, tau, data[2 + 5 * deg : 2 + 6 * deg]),
            )

    def record(self, et, target, observer):
        # TODO: improve search algorithm if this is too slow
        for (etbeg, etend, t, o, frame, type, rbeg, rend)  in self._summary_records:
            if (t, o) == (target, observer):
                if (etbeg <= et <= etend):
                    return self.fetch_record(et, type, rbeg, rend)
        raise RuntimeError("No record found covering that et/target/observer combination")

    def fetch_record(self, et, type, rbeg, rend):
        offset = (rend - 4) * 8 # the last four words
        fmt = "dddd"
        init, intlen, rsize, n = struct.unpack_from(fmt, self._mem, offset=offset)
        internal_offset = math.floor((et - init) / intlen) * rsize
        record = 8 * int(rbeg + internal_offset)
        if type == 2:
            order = (int(rsize) - 2) / 3 - 1
        elif type == 3:
            order = (int(rsize) - 2) / 6 - 1
        else:
            raise RuntimeError("Only Types I and II are implemented")
        return type, order, array.array("d", self._mem[record - 8 : record + int(rsize) * 8])


if __name__ == "__main__":
    import sys

    if len(sys.argv) != 3:
        print "usage: {0} /path/to/jup310.bsp /path/to/libcspice.so"
        exit()

    spk = SPK(sys.argv[1])

    # generate states of all Galilean Satellites for every hour in year 2000
    ts = [3600.0 * h for h in xrange(366 * 24)]

    from timeit import default_timer as timer
    tbeg = timer()
    ss = [[spk.state(t, s, 5) for s in [501, 502, 503, 504]] for t in ts]
    tend = timer()
    n = 366 * 24 * 4
    print "our performance: {0:,.0f} six-state interpolations per second".format(n / (tend - tbeg))

    # let's test whether we are correct... using SPICE
    import ctypes # Python Standard Library for calling C-libraries
    lib = ctypes.CDLL(sys.argv[2])
    lib.furnsh_c(ctypes.c_char_p(sys.argv[1]))
    lib.trcoff_c() # makes CSPICE faster by turning tracing off
    obs = ctypes.c_int(5)
    def state(t, s):
        out =  (ctypes.c_double * 6)()
        t = ctypes.c_double(t)
        ref = ctypes.c_char_p("J2000")
        lt = ctypes.c_double()
        lib.spkgeo_c(s, t, "J2000", 5, out, ctypes.byref(lt))

        return [i for i in out]
    tbeg = timer()
    sp = [[state(t, s) for s in [501, 502, 503, 504]]for t in ts]
    tend = timer()
    print "cspice performance: {0:,.0f} six-state interpolations per second".format(n / (tend - tbeg))
    def diff(s1, s2):
        return math.sqrt(sum((a - b)**2 for (a, b) in zip(s1, s2)))

    max_diff = 0.0
    for row in xrange(366 * 24):
        for col in xrange(4):
            d = diff(ss[row][col], sp[row][col])
            if (d > max_diff):
                a = ss[row][col]
                b = sp[row][col]
                max_diff = d

    print "Max diff: {0:23.16e}\n".format(max_diff)
    print a
    print b

Sample run

python spk.py /data/spice/jup310.bsp /usr/local/lib/libcspice.so 
our performance: 73,604 six-state interpolations per second
cspice performance: 238,672 six-state interpolations per second
Max diff:  3.4924596548334554e-10

(1232518.4345062333, -1281302.1300380777, -587586.9117965015, 6.1812059185567, 4.901470143025447, 2.4082200089953774)
[1232518.434506233, -1281302.130038078, -587586.9117965016, 6.181205918556701, 4.901470143025446, 2.408220008995377]

Hi! I would like to convert the jup310.bsp to the "usual" (ascii or binary) DE planetary ephemeris format (ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/) to read it using the JPL provided routines. What you did here and here (https://space.stackexchange.com/questions/12506/what-is-the-exact-format-of-the-jpl-ephemeris-files) looks like the closest I could find to a solution to my issue. I cannot use SPICE within my soft, I need to import the ChebyCoeff and interpolate them (possibly in the correct way) to get Europa ephemeris.
Thanks for any links on how to perform this format conversion! ^^