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
