Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Example demonstrating how could one work with SPICE SPK files in their original DAF format.
"""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 <Juan.Arrieta@nablazerolabs.com>
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
@arrieta

This comment has been minimized.

Copy link
Owner Author

commented Mar 19, 2017

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]
@steo85it

This comment has been minimized.

Copy link

commented Jun 12, 2017

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! ^^

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.