10
$\begingroup$

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?

$\endgroup$
12
$\begingroup$

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.

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

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 loaded segment
    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)) {
                err("bad segment");
                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) {
                err("target %d not found!", s[i].target);
                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');
        }
        unload_c(*argv);
        spk_free(s, n);
    }
    return 0;
}
$\endgroup$
  • $\begingroup$ This is exactly what I was looking for. Thanks so much, Mark! $\endgroup$ – Doctor Mohawk Aug 6 '15 at 15:57
  • 2
    $\begingroup$ So, having done a bit more research, the SPICE library is too large to include and too tangled to pick out pieces I'd like and use those. Any thoughts? $\endgroup$ – Doctor Mohawk Aug 11 '15 at 19:49
  • 1
    $\begingroup$ @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. $\endgroup$ – Mark Adler Aug 12 '15 at 8:44
  • $\begingroup$ @MarkAdler - Very late, but Doctor Mohawk was wanting something that could be embedded in flight software. The SPICE library is not qualified for that kind of use. It isn't even qualified as NASA class C software, let alone the A or B classification needed for flight software. $\endgroup$ – David Hammen Apr 28 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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