# Is there a way to extract the Chebyshev coefficients for a body from a SPICE kernel?

I want to extract the Chebyshev coefficients for the positions of the Moon and Sun relative to the Earth from a SPICE kernel (mainly because DE430 is too large for me to include in a software package, and I don't need all the other ephemerides, but also because I don't want to include SPICE code in my software release directly).

Any thoughts on what SPICE method I might use to produce these coefficients? My mind is a little boggled by all of the information in these packages.

Or is there a better way that I'm not considering currently?

I want to extract the Chebyshev coefficients for the positions of the Moon and Sun relative to the Earth from a SPICE kernel (mainly because DE430 is too large ...)

You've described almost half (44.8%) of DE430 with just those three items, the positions of the Moon and Sun relative to the Earth. This is not why DE430 is so large. The primary reason DE430 is so large is the large time it covers, 1550 AD to 2650 AD.

With a huge amount of programming work, you can cut the size down by a factor of eleven if you're only interested in 1950 AD to 2050 AD. With a bit more programming work you can cut the size by another factor of two by eliminating everything but the position of the Earth-Moon barycenter, the Moon, and the Sun.

If you go to ftp://ssd.jpl.nasa.gov/pub//eph/planets/ascii/de430 , you'll see a bunch of ASCII files. There are eleven files named ascp[yyyy].430 and two files named header.430_[xxx]. Suppose you're only interested in the period 1950 AD to 2050 AD. You need the Chebychev coefficients in ascp1950.430, plus the stuff in header.430_229 that helps you read those coefficients. If you look at that header file, you'll see that it starts with two coefficients KSIZE= 2036 NCOEFF= 1018 followed by data organized in groups. That NCOEFF= 1018 is important. It indicates the number of Chebychev coefficients per block (more on this later). With regard to the GROUPs,

• GROUP 1030 contains

• the start time (as a TDB Julian date),
• the end time (as a TDB Julian date), and
• the number of days in a block.
• In the case of DE430, it starts on JD 2287184.5, ends on JD 2688976.5, and is organized in blocks of 32 days each.
• GROUP 1040 contains the names of 229 parameters.

• GROUP 1041 contains the values of those 229 parameters. You'll want the value of EMRAT, the ratio of the Earth's mass to that of the Moon. I'll leave it up to you to determine what to do with that, but you truly do want it.

• GROUP 1050 contains 3 lines of 13 integers each.

• The first line contains the initial word (Fortran numbering) of the Chebychev coefficients in a block of coefficients (more on this later) for

1. Mercury,
2. Venus,
3. the Earth-Moon barycenter,
4. Mars,
5. Jupiter,
6. Saturn,
7. Uranus,
8. Neptune,
9. Pluto,
10. the Moon,
11. the Sun,
12. the librations of the Earth, and
13. the mutations of the Moon.
• The second line contains the number of Chebychev polynomial coefficients (one plus the polynomial order) for each those objects.

• The third line contains the number of sets of Chebychev polynomial coefficients per 32-day block for each of those objects.

For example, the parameters for the Earth-Moon barycenter are:

• Initial word = 231
• Number of Chebychev coefficients = 13
• Number of sets of coefficients per block = 2

Every item except for #12 (Earth librations) represents a set of three values. In that special case of Earth librations, there are only two values.

If you look at (for example) ascp1950.430, you'll see that it starts with the line " 1 1018" followed by 1017 numbers. There's a missing value; there should be 1018 numbers. (Oops! Somebody at JPL didn't test their software. That's just one of a number of wonderful features.) You'll need to fix that missing value bug (just add a 0.0 to complete the set). Then you'll need to interpret the values per the GROUP 1050 spec.

For example, the Earth-Moon barycenter comprises two sets of 12th order Chebychev polynomial coefficients for each of x, y, and z, starting at word 231. The first thirteen values are the values for the first set of coefficients for x, and so on.

The first two words in each block of 1018 values are very important. These are the start and end times (in TDB Julian days) of the block. In general, you never, ever want to extrapolate when dealing with Chebychev polynomials. The DE4xx comprises coefficients for sets of piecewise continuous Chebychev polynomials. You absolutely need to find the right set. Those first two words of each block, coupled with the object-specific number of sets of polynomial coefficients per block, will target you toward the correct set of coefficients for the point in time and object of interest.

To use these coefficients, you need to

• Transform your time (presumably something like UTC, TAI, or TT) to TDB.
• Find the block (and the sub-block in the case of objects such as the Earth-Moon barycenter) that contains the coefficients pertinent to that object and that point in time,
• Retrieve the Chebychev coefficients for that time interval,
• Offset and scale the point in time to a value between -1 (start of interval) to +1 (end of interval), and finally,
• Compute the values of the x,y, and z Chebychev polynomials given that scaled/offset time value.

Alternatively, you could just use the SPICE libraries and SPICE kernels. I went through this mess over a decade ago. The effort to do so represents a few months of my life that I want back.

• Dave, do you ever sleep? Thanks for this. I just want to give you all of my karma after reading it. – Autumnsault Aug 6 '15 at 14:12
• Amazing. Thanks! This also seems helpful, providing much more detailed explanations and IDL code: asd.gsfc.nasa.gov/Craig.Markwardt/bary – nealmcb Jun 4 '17 at 20:08

What's wrong with including some of the SPICE code? It is a library, so it would only link in what you need. SPICE can be used freely in any products, including commercial, for-profit use.

Ideally you would still use the SPICE code to get states from the kernel. Then you can just edit down the kernel to include only the bodies and time spans of interest. SPKMERGE will do that for you. With this input to SPKMERGE:

LEAPSECONDS_KERNEL = naif0010.tls
SPK_KERNEL = pared-down.bsp
SOURCE_SPK_KERNEL = de430.bsp
BODIES = 10 301 399
BEGIN_TIME = 2001-01-01
END_TIME = 2051-01-01


the resulting pared-down.bsp with only the Moon, Earth, and the Sun over the specified date range, is 3.2 MB, as compared to de430.bsp which is 114 MB.

Here is a program that extracts the Chebyshev coefficients from an SPK file into a memory structure that is then used to compute positions and velocities without use of the SPICE library.

/* cheby.c - extract and use Chebyshev position-only coefficients from an SPK
file.  This illustrates how the coefficients may be extracted once using the
SPICE library, and then used after that without the SPICE library.

Mark Adler   August 15, 2015   placed in the public domain
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <math.h>
#include "SpiceUsr.h"

/*
DAF/SPK file format notes:

ic[0] - target code
ic[1] - center code
ic[2] - frame code
ic[3] - representation code (2 == Chebyshev position only)
ic[4] - initial address of array
ic[5] - final address of array

len = ic[5] - ic[4] + 1

dc[0] - initial epoch of data (seconds relative to J2000)
dc[1] - final epoch of data (seconds relative to J2000)

seg[len-4] - initial epoch of the first record (seconds relative to J2000)
seg[len-3] - interval length of each record (seconds)
seg[len-2] - elements in each record
seg[len-1] - number of records

seg[len-1] * seg[len-2] + 4 == len
n = seg[len-2]
num = (n - 2) / 3

rec[0] - midpoint of interval covered by record (seconds relative to J2000)
rec[1] - radius of interval (seconds)
rec[2..num+1] - X coefficients (constant term first)
rec[num+2..2*num+1] - Y coefficients
rec[2*num+2..n-1] - Z coefficients

For t, evaluate the Chebyshev polynomials T_n at (t - rec[0]) / rec[1],
multiply by the coefficients, and sum.  The derivatives of the polynomials
can be used to compute the velocity.  See cheby_eval() here.  The results
are in km and km/s.  Note that all times are ephemeris times, and so do not
take into account leap seconds.

A SpiceDouble is simply a C double.  A SpiceInt is an integer type whose
size is half that of double, so that two SpiceInt's fit in a SpiceDouble.
*/

/* Evaluate the given Chebyshev polynomial at x, returning both the evaluated
polynomial in *f, and the evaluated derivative of the polymonial in *df. The
number of coefficients is num (the degree of the polynomial is num - 1), and
the coefficients are coeff[0..num-1].  The first coefficient coeff[0] is the
constant term.  The scaling of x is provided by the midpoint scale[0] and
the radius scale[1].  x must fall in the range scale[0] - scale[1] to
scale[0] + scale[1].  Outside of that range, the polynomial is not valid. */
void cheby_eval(double x, double *scale, double *coeff, long num,
double *f, double *df)
{
double x2, w0 = 0., w1 = 0., dw0 = 0., dw1 = 0., tmp;

x = (x - scale[0]) / scale[1];
x2 = x * 2.;
while (--num) {
tmp = dw1;
dw1 = dw0;
dw0 = w0 * 2. + dw0 * x2 - tmp;
tmp = w1;
w1 = w0;
w0 = coeff[num] + (x2 * w0 - tmp);
}
*f = coeff[0] + (x * w0 - w1);
*df = (w0 + x * dw0 - dw1) / scale[1];
}

/* Find the appropriate SPK record for time t and compute the position and
velocity for that time.  Returns 0 on success, 1 if the time is not covered
by the segment. */
int cheby_posvel(double t, double *seg, long len, double pos[3], double vel[3])
{
long k, num;

k = (long)floor((t - seg[len - 4]) /    // seg[len-4] is initial epoch
seg[len - 3]);          // seg[len-3] is record span
if (k < 0 || k >= (long)seg[len - 1])   // seg[len-1] is number of records
return 1;
num = (long)seg[len - 2];               // seg[len-2] is size of record
seg += k * num;                         // point seg to the record for t
num = (num - 2) / 3;                    // number of coefficients
cheby_eval(t, seg, seg + 2, num, pos, vel);
cheby_eval(t, seg, seg + 2 + num, num, pos + 1, vel + 1);
cheby_eval(t, seg, seg + 2 + 2 * num, num, pos + 2, vel + 2);
return 0;
}

/* Verify that the provided segment meets the constraints of a uniform set of
coefficient records.  Return 0 on success or 1 if the segment is invalid.
This should be done before using the segment in order to avoid segfaults on
invalid data. */
int cheby_verify(double *seg, long len)
{
double
recs = seg[len - 1],    // number of records
elts = seg[len - 2],    // elements (doubles) in each record
span = seg[len - 3],    // time span of each record in seconds
init = seg[len - 4];    // initial epoch in seconds relative to J2000
long n, k;
double *p, *q;

if (recs != (long)recs ||                       // recs is an integer
elts != (long)elts ||                       // elts is an integer
(long)recs * (long)elts + 4 != len ||       // total length is correct
3 * (((long)elts - 2) / 3) + 2 != elts ||   // integer number of coeffs
seg[0] - seg[1] != init ||                  // 1st start is init
span != 2 * seg[1])                         // 1st radius matches span
return 1;
n = (long)recs;
k = (long)elts;
p = seg;
while (--n) {
q = p + k;                                  // scan all q following p
if (q[1] != p[1] ||                         // all radii the same
q[0] - q[1] != p[0] + p[1])             // next start is last end
return 1;
p = q;
}
return 0;
}

/* Print an error message. */
void err(char *msg, ...)
{
fputs("cheby error: ", stderr);
va_list ap;
va_start(ap, msg);
vfprintf(stderr, msg, ap);
va_end(ap);
putc('\n', stderr);
}

/* SPK segment descriptor. */
typedef struct {
long target;        // target body code
long center;        // center body code
long frame;         // frame of reference code
long len;           // length of segment in doubles
double *seg;        // allocated segment
} segment_t;

/* Load one segment of an SPK file, which covers one target over a range of
epochs.  Save the target code, reference location code for the target
position, and the reference frame code.  Load the segment and verify its
structure.  On success return 0.  If there is an error, return 1 and set
s->seg to NULL. */
int cheby_segment(SpiceInt daf, SpiceDouble *dc, SpiceInt *ic, segment_t *s)
{
SpiceDouble *last;

// save segment codes
s->target = ic[0];
s->center = ic[1];
s->frame = ic[2];

// allocate memory for the segment and read it in
s->len = ic[5] - ic[4] + 1;                 // number of doubles in segment
s->seg = malloc(s->len * sizeof(SpiceDouble));
if (s->seg == NULL) {
err("out of memory");
return 1;
}
dafgda_c(daf, ic[4], ic[5], s->seg);        // load segment
if (failed_c()) {
reset_c();
free(s->seg);
s->seg = NULL;
err("could not read SPK segment from file");
return 1;
}

// verify the integrity of the segment
last = s->seg + s->len - 4 - (long)(s->seg[s->len - 2]);
if (cheby_verify(s->seg, s->len) ||         // segment structure ok
dc[0] != s->seg[s->len - 4] ||          // start epoch matches
dc[1] != last[0] + last[1]) {           // end epoch matches
free(s->seg);
s->seg = NULL;
err("SPK segment format is invalid");
return 1;
}

return 0;
}

/* Scan through the SPK file path and extract all of the Chebyshev
position-only segments, saving them in an allocated array of segment_t,
which is returned.  If there is an error, NULL is returned.  *segs is set to
the number of segments in the array.  Once this is done, this array can be
used by cheby_verify() and cheby_posvel() above, with no dependency on or
reference to the SPICE library. */
segment_t *spk_extract(char *path, long *segs)
{
SpiceInt daf;
SpiceBoolean found;
union {
SpiceDouble d[128];
SpiceChar c[1024];
} sum;
const SpiceInt nd = 2, ni = 6;
SpiceDouble dc[nd];
SpiceInt ic[ni];
segment_t *spk, *mem;

// turn off error reporting and aborts for SPICE functions
errprt_c("set", 0, "none");
erract_c("set", 0, "return");

// open the file and verifiy that it is a DAF SPK file
dafopr_c(path, &daf);                       // open SPK file for reading
if (failed_c()) {
reset_c();
err("could not open %s as a DAF", path);
return NULL;
}
dafgsr_c(daf, 1, 1, 128, sum.d, &found);    // read first record
if (failed_c() || !found || memcmp(sum.c, "DAF/SPK ", 8)) {
reset_c();
dafcls_c(daf);
err("%s is not an SPK file", path);
return NULL;
}

// count the number of Chebyshev position-only segments in the DAF file
*segs = 0;
dafbfs_c(daf);                              // begin forward search
while (daffna_c(&found), found) {           // find the next array
dafgs_c(sum.d);                         // get array summary
dafus_c(sum.d, nd, ni, dc, ic);         // unpack the array summary
if (failed_c())
break;
if (ic[3] == 2)                         // Chebyshev position only
(*segs)++;                          // count segment
}
if (failed_c() || *segs == 0) {
reset_c();
dafcls_c(daf);
err("file error or Chebyshev position-only segments in %s", path);
return NULL;
}

// allocate table of segment descriptors
spk = malloc(*segs * sizeof(segment_t));
if (spk == NULL) {
dafcls_c(daf);
err("out of memory");
return NULL;
}

// read and save the Chebyshev position-only segments
*segs = 0;
dafbfs_c(daf);                              // begin forward search
while (daffna_c(&found), found) {           // find the next array
dafgs_c(sum.d);                         // get array summary
dafus_c(sum.d, nd, ni, dc, ic);         // unpack the array summary
if (failed_c())
break;
if (ic[3] == 2 && !cheby_segment(daf, dc, ic, spk + *segs))
(*segs)++;
}
if (failed_c() || *segs == 0) {
reset_c();
dafcls_c(daf);
free(segs);
err("no valid Chebyshev position-only segments in %s", path);
return NULL;
}

// close the DAF file and return segment table
dafcls_c(daf);
errprt_c("set", 0, "short");
erract_c("set", 0, "abort");
mem = realloc(spk, *segs * sizeof(segment_t));
if (mem != NULL)
spk = mem;
return spk;
}

/* Free the resources of an SPK structure created by spk_extract(). */
void spk_free(segment_t *s, long n)
{
long i;

for (i = 0; i < n; i++)
free(s[i].seg);
free(s);
}

/* Load the SPK files on the command line and verify the position and velocity
at J2000 + 0 seconds for each Chebyshev position-only segment against the
same result from the SPICE library. */
int main(int argc, char **argv)
{
segment_t *s;
long n, i;
double pos[3], vel[3];
SpiceInt eph, frame, center;
SpiceDouble desc[5], pv[6];
SpiceBoolean found;
SpiceChar id[41];

while (++argv, --argc) {
s = spk_extract(*argv, &n);
if (s == NULL) {
err("could not load %s as an SPK file", *argv);
continue;
}
furnsh_c(*argv);
for (i = 0; i < n; i++) {
// show segment info and position and velocity at J2000 + 0
printf("target = %ld, center = %ld, frame = %ld\n",
s[i].target, s[i].center, s[i].frame);
if (s[i].seg == NULL || cheby_verify(s->seg, s->len)) {
putchar('\n');
continue;
}
if (cheby_posvel(0, s[i].seg, s[i].len, pos, vel)) {
err("J2000 + 0 out of range (!)");
putchar('\n');
continue;
}
printf("pos(0) = (%g, %g, %g)\n", pos[0], pos[1], pos[2]);
printf("vel(0) = (%g, %g, %g)\n", vel[0], vel[1], vel[2]);

// check position and velocity against SPICE library access
spksfs_c(s[i].target, 0, sizeof(id), &eph, desc, id, &found);
if (!found) {
putchar('\n');
continue;
}
spkpvn_c(eph, desc, 0, &frame, pv, &center);
if (s[i].frame != frame || s[i].center != center)
err("codes mismatch");
if (pos[0] != pv[0] || pos[1] != pv[1] || pos[2] != pv[2])
err("position mismatch");
if (vel[0] != pv[3] || vel[1] != pv[4] || vel[2] != pv[5])
err("velocity mismatch");
putchar('\n');
}

• @Dr.JohnnyMohawk: if you must, I have updated the answer with a coefficient extractor. The very simple chbint() routine in the SPICE source code uses the coefficients to get a position component and a velocity component. – Mark Adler Aug 12 '15 at 8:44