ephcom.c

//   so the JPL software uses that value because it reads in Fortran 'A6'
//   character strings.
//
    fprintf( outfp, "%-84s%-84s%-84s", header->ttl[0], header->ttl[1], header->ttl[2] );
    blockout = 3 * EPHCOM_MAXTTL; // Just wrote 3 84-byte strings to start output file
//
//   Now output 400 cnam entries to the output file.
//
    for ( i = 0; i < header->ncon; i++ )
    {
        fprintf( outfp, "%-6s", header->cnam[i] );
        blockout += 6;
    }
    for (; i < 400; i++ )
    {
        fprintf( outfp, "      " ); // Round out to 400 entries, all blank at end
        blockout += 6;
    }
//
//   Binary values: Make sure bytes are in big-endian (network) order for file.
//
    for ( i = 0; i < 3; i++ )
    {
        ephcom_outdouble( outfp, header->ss[i] ); // Write net-order bytes from double precision
        blockout += 8;
    }
    ephcom_outint( outfp, header->ncon );
    blockout += 4;
    ephcom_outdouble( outfp, header->au );
    blockout += 8;
    ephcom_outdouble( outfp, header->emrat );
    blockout += 8;
//
//   If there are no coefficients for an ipt[i][] object (i.e., ipt[i][1]==0),
//   then ipt[i][0] should contain the value of the next available coefficient
//   number rather than 0, as per communication of Myles Standish to Paul Hardy
//   on preferred format of ephemeris headers.
//
//   If there are no libration coefficients (i.e., lpt[1]==0), then lpt[0]
//   should contain the value of the next available coefficient number rather
//   than 0 as well, as per the same communication from Myles Standish.
//
// First set j to maximum index into ipt[] that has coefficients
    j = 0;
    for ( i = 1; i < 12; i++ )
        if ( header->ipt[i][1] > 0 && header->ipt[i][0] > j )
            j = i;
// Now set j to next available index count.
    if ( header->lpt[1] > 0 && header->lpt[0] > j )
        j = header->lpt[1] + header->lpt[1] * header->lpt[2] * 3;
    else
        j = header->ipt[j][0] +
            header->ipt[j][1] * header->ipt[j][2] * ( j == 11 ? 2 : 3 );
    for ( i = 1; i < 12; i++ )
        if ( header->ipt[i][0] == 0 )
            header->ipt[i][0] = j;
    if ( header->lpt[0] == 0 )
        header->lpt[0] = j;

    for ( j = 0; j < 12; j++ )
    {
        for ( i = 0; i < 3; i++ )
        {
            ephcom_outint( outfp, header->ipt[j][i] );
            blockout += 4;
        }
    }
    ephcom_outint( outfp, header->numde );
    blockout += 4;
    for ( i = 0; i < 3; i++ )
    {
        ephcom_outint( outfp, header->lpt[i] );
        blockout += 4;
    }
//
//   Now pad the end of the first record with null bytes.  Note: the
//   JPL Fortran software just skips to next record at this point.
//
    for ( i = blockout; i < blockbytes; i++ )
    {
        fputc( '\0', outfp );
    }
//
//   End of first block.  Now set blockout to 0 and start with next block.
//
    blockout = 0;
    for ( i = 0; i < header->ncon; i++ )
    {
        ephcom_outdouble( outfp, header->cval[i] );
        blockout += 8;
    }
//
//   Pad with double-precision zeroes for rest of array.
//
    for (; i < 400; i++ )
    {
        ephcom_outdouble( outfp, (double) 0.0 );
        blockout += 8;
    }
//
//   Pad with nulls for rest of block.
//
    for ( i = blockout; i < blockbytes; i++ )
    {
        fputc( '\0', outfp );
    }
//
//   Finished normally.
//
}

//! Write JPL ephemeris coefficient block of data in binary form.
//! @param outfp [IN ONLY]Pointer to the FILE which will be used to
//! output the binary block of data.
//! @param header [IN ONLY]Pointer to an ephcom_Header struct which
//! contains the JPL ephemeris header information.  Some of these
//! data are used to control where in the binary FILE the block of
//! data is written and how large it is.
//! @param blocknum [IN ONLY]Block number.  This value helps control
//! where in the binary FILE the block of data is written.
//! @param datablock [IN ONLY]Pointer to an array that contains the
//! block of data that will be output in binary form.
//!
void
ephcom_writebinary_block(
    FILE * outfp,
    ephcom_Header *header,
    int blocknum,
    double *datablock )
{
    int i;
    int filebyte;
    int filepos;

//
//   Find out where we need to point in the binary file.
//
    filebyte = ( blocknum + 2 ) * header->ncoeff * 8; // 8 bytes per coefficient
//
//   If the file isn't that large, pad it with null bytes
//
    fseek( outfp, 0L, SEEK_END );
    filepos = ftell( outfp );
    if ( filepos < filebyte )
    {
        for ( i = 0; i < ( filebyte - filepos ); i++ )
        {
            fputc( '\0', outfp );
        }
    }
//
//   Now go to position where we want to start writing.
//
    fseek( outfp, filebyte, SEEK_SET );
    for ( i = 0; i < header->ncoeff; i++ )
    {
        ephcom_outdouble( outfp, datablock[i] );
    }
}


//! Write out a block of JPL binary data in ascii form that is nicely
//! formatted and therefore verbose.
//!
//! @param outfp [IN ONLY]Pointer to the FILE which will be used to
//! output the nicely formatted data.
//! @param header [IN ONLY]Pointer to an ephcom_Header struct which
//! contains the JPL ephemeris header information.  Some of these
//! data are used to control how to output the block of data in a
//! nicely formatted way.
//! @param datablock [IN ONLY]Pointer to an array that contains the
//! block of data that will be output in nicely formatted form.
//!
void
ephcom_parse_block(
    FILE * outfp,
    ephcom_Header *header,
    double *datablock )
{
    int i0, i1, i2, i3;
    int blockword;
//
//   Names of the objects in Chebyshev coefficient arrays.
//
    const char *ephcom_coeffname[13] = {
        "Mercury", "Venus",           "EMBary", "Mars",     "Jupiter", "Saturn", "Uranus", "Neptune",
        "Pluto",   "Geocentric Moon", "Sun",    "Nutation", "Libration"
    };

    blockword = 0;
    fprintf( outfp, "@%04d StartJD\t%12.2f\n", blockword++, datablock[0] );
    fprintf( outfp, "@%04d EndJD\t%12.2f\n", blockword++, datablock[1] );
    for ( i0 = 0; i0 < 13; i0++ ) // For all bodies
    {
        fprintf( outfp, "Body\t%d (%s)\n", i0 + 1, ephcom_coeffname[i0] );
        for ( i1 = 0; i1 < header->ipt[i0][2]; i1++ ) // For all subintervals
        {
            fprintf( outfp, "  Subinterval %d of %d\n", i1 + 1, header->ipt[i0][2] );
            for ( i2 = 0; i2 < ( i0 == 11 ? 2 : 3 ); i2++ ) // For all coordinates
            {
                fprintf( outfp, "    %cCoefficients\n", 'X' + i2 );
                for ( i3 = 0; i3 < header->ipt[i0][1]; i3++ ) // For all coefficients
                {
                    blockword = header->ipt[i0][0] +
                                i1 * header->ipt[i0][1] * ( i0 == 11 ? 2 : 3 ) +
                                i2 * header->ipt[i0][1] + i3 - 1;
                    fprintf( outfp, "      @%04d [%2d of %2d] %25.17E\n",
                        blockword, i3 + 1, header->ipt[i0][1], datablock[blockword] );
                }
            }
        }
    }
}

//! Read the next two lines of the ascii form of a JPL ephemeris file.
//! For the last of those check that it is exactly the same as the
//! expected group header form, that is, 12 characters containing
//! "GROUP nnnn" with the same nnnn as expected.  If the group is not
//! what is expected, then print an error message and exit.
//!
//! @param group [OUT ONLY]Pointer to a group header that upon return
//! will be filled by a 12-character null-terminated string that we
//! read.
//! @param expected [IN ONLY]Pointer to a 12-character null-terminated
//! string that we expect to read.
//! @param infile [IN ONLY]Pointer to the JPL ascii ephemeris FILE
//! that we read.
//!
static void
ephcom_nxtgrp( char *group, const char *expected, FILE *infile )
{
    char readbuf[EPHCOM_MAXLINE + 1];

    fgets( readbuf, EPHCOM_MAXLINE, infile ); // Blank Line
    fgets( readbuf, EPHCOM_MAXLINE, infile ); // "GROUP   dddd\n"
    strncpy( group, readbuf, 12 );
    group[12] = '\0';
    if ( strncmp( group, expected, 12 ) != 0 )
    {
        fprintf( stderr, "Badly formed header; \"%s\" not found.\n\n", expected );
        exit( 1 );
    }
    fgets( readbuf, EPHCOM_MAXLINE, infile ); // Blank Line
}

//! Write a double-precision value to the given binary file with
//! bytes swapped if necessary to match network order (Big Endian).
//! On Intel 80x86 the bytes will get swapped, on Motorola or SPARC
//! they won't.
//!
//! @param outfp [IN ONLY]Pointer to the binary output FILE.
//! @param x [IN ONLY]Double-precision value that will be written to
//! the binary FILE.
//!
static void
ephcom_outdouble( FILE *outfp, double x )
{
    double        retval;
    unsigned char ch[8];

    memcpy( (void *) ch, (const void *) &x, 8 );
    gnulliver64c( ch );
    fwrite( ch, 1, 8, outfp );
}

//! Write a integer value to the given binary file with bytes swapped
//! if necessary to match network order (Big Endian).  On Intel 80x86
//! the bytes will get swapped, on Motorola or SPARC they won't.
//!
//! @param outfp [IN ONLY]Pointer to the binary output FILE.
//! @param u [IN ONLY]Integer value that will be written to the binary
//! FILE.
//!
static void
ephcom_outint( FILE * outfp, unsigned u )
{
    unsigned u2;

    u2 = gnulliver32( u );
    fwrite( &u2, 4, 1, outfp );
}

//! Read a double-precision value from the given binary file with
//! bytes swapped if necessary to match host endian order.  On Intel
//! 80x86 the bytes will get swapped, on Motorola or SPARC they won't.
//!
//! @param infp [IN ONLY]Pointer to the binary input FILE.
//! @returns the (possibly) byte-swapped double-precision value that
//! was read from the binary input FILE.
//!
double
ephcom_indouble( FILE *infp )
{
    double        x;
    double        retval;
    unsigned char ch[8];

    //
    // Handle as character string until bytes are in correct order,
    // then copy to double once they are.
    //
    fread( ch, 1, 8, infp );
    gnulliver64c( ch );
    memcpy( (void *) &retval, (const void *) ch, (size_t) 8 );
    return ( retval );
}

//! Read an integer value from the given binary file with bytes
//! swapped if necessary to match host endian order.  On Intel 80x86
//! the bytes will get swapped, on Motorola or SPARC they won't.
//!
//! @param infp [IN ONLY]Pointer to the binary input FILE.
//! @returns the (possibly) byte-swapped integer value that was read
//! from the binary input FILE.
//!
int
ephcom_inint( FILE *infp )
{
    unsigned u;
    int      retval;

    fread( &u, 4, 1, infp );
    retval = (int) gnulliver32( u );
    return ( retval );
}

//! Function to convert a string with a double precision value written
//! in C to a double precision value that fortran understands (i.e.,
//! with a "D" exponent character).  Conversion happens in place.
//!
//! @param buf [IN AND OUT]Pointer to a character string holding the C
//! double-precision value (in "E" exponential format with a trailing
//! blank after the exponent to make room for the leading zero
//! notation of the fortran value) on input and the fortran
//! double-precision value in "D" exponential format for fortran (with
//! leading zero before the decimal point) on output.  If the resulting
//! fortran exponent has a leading zero and more than two digits that
//! leading zero is dropped while right justification is maintained
//! by shifting the whole string to the right by one byte to overlay
//! that leading zero by the exponent sign.
//!
static void
ephcom_doublestrc2f( char *buf )
{
    int    i, j, istart, istop, exp, edigits;
    double x;

    // Deal with three or more digit exponent with leading zero
    // (which can occur for the Windows case).
    for ( istop = 0; toupper( buf[istop] ) != 'E'; istop++ )
        ;

    // buf[istop] is 'E', buf[istop+1] is the exponent sign, buf[istop+2],...
    // is the absolute value of the exponent.
    if ( buf[istop + 2] == '0' && isdigit( buf[istop + 3] ) && isdigit( buf[istop + 4] ) )
    {
        for ( istart = istop + 2; istart > 0; istart-- )
            buf[istart] = buf[istart - 1];
        buf[0] = ' ';
        istop++;
    }

    for ( istart = 0; isspace( buf[istart] ); istart++ )
        ;
    x = atof( &buf[istart] );

    exp = atoi( &buf[istop + 1] );
    exp++;
    if ( exp < 0 )
    {
        buf[istop + 2] = '-';
        exp            = -exp;
    }
    else
    {
        buf[istop + 2] = '+';
    }
    if ( x == 0.0 )
        exp = 0;
    if ( exp < 100 )
        edigits = 2;
    else if ( exp < 1000 )
        edigits = 3;
    else
        edigits = 4;

    while ( edigits > 0 )
    {
        buf[istop + edigits + 2] = exp % 10 + '0';
        exp /= 10;
        edigits--;
    }

    buf[istop + 1] = 'D';

    while ( istop > istart && buf[istop - 1] != '.' )
    {
        buf[istop] = buf[istop - 1];
        istop--;
    }

    buf[istop]     = buf[istop - 2]; // buf[istop-1] == '.'
    buf[istop - 2] = '0';            // leading zero
}

//! This ephcom_pleph routine takes coordinates already calculated in
//! a ephcom_Coords struct and returns selected results in an array
//! depending on the ntarg and ncentr indices provided by the calling
//! routine.  These indices have the following interpretation (note
//! the offset of one compared to the interpretation of the first index
//! of pv in the ephcom_Coords struct.
//! <table border>
//!   <tr> <td><b>Index</b></td> <td><b>Identification</b></td> </tr>
//!   <tr> <td>1</td> <td>Mercury</td> </tr>
//!   <tr> <td>2</td> <td>Venus</td> </tr>
//!   <tr> <td>3</td> <td>Earth</td> </tr>
//!   <tr> <td>4</td> <td>Mars</td> </tr>
//!   <tr> <td>5</td> <td>Jupiter</td> </tr>
//!   <tr> <td>6</td> <td>Saturn</td> </tr>
//!   <tr> <td>7</td> <td>Uranus</td> </tr>
//!   <tr> <td>8</td> <td>Neptune</td> </tr>
//!   <tr> <td>9</td> <td>Pluto</td> </tr>
//!   <tr> <td>10</td> <td>Moon</td> </tr>
//!   <tr> <td>11</td> <td>Sun</td> </tr>
//!   <tr> <td>12</td> <td>Solar System Barycenter</td> </tr>
//!   <tr> <td>13</td> <td>Earth-Moon Barycenter</td> </tr>
//!   <tr> <td>14</td> <td>Nutation Angles</td> </tr>
//!   <tr> <td>15</td> <td>Libration Angles</td> </tr>
//!   <tr> <td>16</td> <td>Moon (Geocentric)</td> </tr>
//! </table>
//!
//! @param coords [IN ONLY]Pointer to a ephcom_Coords struct which
//! contains interpolated values of all coordinates and their time derivatives
//! as calculated from Chebyshev coefficients supplied by a JPL ephemeris.
//! @param ntarg [IN ONLY]Index interpreted according to the
//! above table which identifies the "target" data in coords.
//! @param ncntr [IN ONLY]Index interpreted according to the above
//! table which identifies the "center" data in coords.  If either
//! ntarg or ncent is 14, the interpolated nutation angle data (two
//! angles and two angle time derivatives) are returned in r.
//! Otherwise, if either ntarg or ncent is 15, the interpolated
//! libration angle data (three angles and three angle time
//! derivatives) are returned in r.  Otherwise, if either ntarg or
//! ncent is 16, the moon geocentric position and velocity data (3
//! positions and 3 velocities) are returned in pv.  Otherwise (the
//! normal case) if 0 < ntarg < 14 and 0 < ncent < 14, the
//! interpolated positions and velocities corresponding to the center
//! index are subtracted from the interpolated positions and
//! velocities corresponding to the target index and the resulting
//! data (3 positions and 3 velocities) are returned in pv.
//! Otherwise, ephcom_pleph() issues an error message and exits.
//! @param r [OUT ONLY]Pointer to an array which upon return will be
//! filled with the requested interpolated results (4 values for
//! nutation and 6 values for everything else).
//!

void
ephcom_pleph( ephcom_Coords *coords, int ntarg, int ncntr, double *r )
{
    int i;
    if ( ntarg == 14 || ncntr == 14 )
    {
        for ( i = 0; i < 4; i++ )
            r[i] = coords->pv[13][i];
    }
    else if ( ntarg == 15 || ncntr == 15 )
    {
        for ( i = 0; i < 6; i++ )
            r[i] = coords->pv[14][i];
    }
    else if ( ntarg == 16 || ncntr == 16 )
    {
        for ( i = 0; i < 6; i++ )
            r[i] = coords->pv[15][i];
    }
    else if ( ( 0 < ntarg && ntarg < 14 ) && ( 0 < ncntr && ncntr < 14 ) )
    {
        for ( i = 0; i < 6; i++ )
            r[i] = coords->pv[ntarg - 1][i] - coords->pv[ncntr - 1][i];
    }
    else
    {
        fprintf( stderr, "ephcom_pleph: Invalid combination of ntarg = %d and ncntr = %d\n", ntarg, ncntr );
        exit( EXIT_FAILURE );
    }
}

//! Interpolate positions and velocities of all stored JPL "bodies" at given time from a JPL
//! binary ephemeris file.  Transform the interpolated data to the following index scheme
//! for JPL bodies where the index is one less than the index used for the arguments
//! of ephcom_pleph():
//! <table border>
//!   <tr> <td><b>Index</b></td> <td><b>Identification</b></td> </tr>
//!   <tr> <td>0</td> <td>Mercury</td> </tr>
//!   <tr> <td>1</td> <td>Venus</td> </tr>
//!   <tr> <td>2</td> <td>Earth</td> </tr>
//!   <tr> <td>3</td> <td>Mars</td> </tr>
//!   <tr> <td>4</td> <td>Jupiter</td> </tr>
//!   <tr> <td>5</td> <td>Saturn</td> </tr>
//!   <tr> <td>6</td> <td>Uranus</td> </tr>
//!   <tr> <td>7</td> <td>Neptune</td> </tr>
//!   <tr> <td>8</td> <td>Pluto</td> </tr>
//!   <tr> <td>9</td> <td>Moon</td> </tr>
//!   <tr> <td>10</td> <td>Sun</td> </tr>
//!   <tr> <td>11</td> <td>Solar System Barycenter</td> </tr>
//!   <tr> <td>12</td> <td>Earth-Moon Barycenter</td> </tr>
//!   <tr> <td>13</td> <td>Nutation Angles</td> </tr>
//!   <tr> <td>14</td> <td>Libration Angles</td> </tr>
//!   <tr> <td>15</td> <td>Moon (Geocentric)</td> </tr>
//! </table>
//!
//! @param infp [IN ONLY]Pointer to a binary version of a JPL
//! ephemeris FILE.
//! @param header [IN ONLY]Pointer to an ephcom_Header struct which
//! contains the JPL ephemeris header information that has already
//! been read from that binary JPL ephemeris FILE.
//! @param coords [IN AND OUT]Pointer to a ephcom_Coords struct which
//! contains on input coords->et2 (the split Julian date where the
//! interpolation should occur) and also the control information
//! coords->km (if nonzero, solar system body coordinates [but not
//! nutation and libration angles which are always returned in
//! radians] will be be returned in kilometers rather than
//! astronomical units); coords->second (if nonzero, all solar system
//! body and nutation and libration angle time derivatives will be
//! returned in per second units rather than in per day units);
//! coords->bary (if nonzero, all solar system body coordinates will
//! be relative to the solar system barycenter, otherwise they will be
//! relative to the solar system body indicated by the coords->center
//! index; and coords->center (the index in the range from 0 to 12
//! which cooresponds to a particular solar system body indicated in
//! the above table, and which only needs to be supplied if
//! coords->bary is zero).  Upon return this struct contains
//! interpolated values of all coordinates and their time derivatives
//! as calculated from the Chebyshev coefficients in datablock using
//! the above body index scheme.
//! @param datablock [OUT ONLY]Pointer to an array that upon return
//! will contain the block of Chebyshev coefficient data that has been
//! read that has a Julian date range that contains the specified
//! Julian date in coords where the interpolation occurred.
//! @returns 0 on success and -1 on failure (due to specified Julian
//! date in coords out of range or some i/o error).
//!
int
ephcom_get_coords( FILE * infp,
                   ephcom_Header *header,
                   ephcom_Coords *coords,
                   double *datablock )
{
    double et2[2], pjd[2]; // Ephemeris times, split into coarse (whole) and fine time  in JD
    double filetime;       // JDs since start of ephemeris file
    double blocktime[2];   // JDs since start of data block
    double subtime;        // JDs since start of subinterval in block

    int    i, j, k;
    int    blocknum;
    // Number of subintervals in data block for this body
    int    nsub;
    // Number of subinterval for this body
    int    subinterval;
    // Offset in datablock for current body and subinterval
    int    dataoffset;
    // Span of one subinterval in days
    double subspan;
    // Span of one subinterval in days (coords->second is zero) or
    // seconds (coords->second is nonzero) used for normalization
    // of the Chebyshev polynomial time derivatives.
    double norm_span;
    // Normalized Chebyshev time, in interval [-1,1].
    double chebytime;
    // Number of coordinates for position and velocity
    int    ncoords;
    // Number of Chebyshev coefficients per coordinate
    int    ncf;
    // Return value
    int    retval;

    // Assume normal return
    retval = 0;

    // et2 is a transformed version of coords->et2 such that the first
    // value of et2 is exactly half-integral to match the exactly
    // half-integral characteristics of the Julian dates returned from
    // the binary ephemeris.
    et2[1] = ephcom_split( coords->et2[0] - 0.5, &et2[0] );
    pjd[1] = ephcom_split( coords->et2[1], &pjd[0] );
    // et2[0] should end up as exactly half integral.
    et2[0] += pjd[0] + 0.5;
    // Deal with fractional remainders.
    et2[1]  = ephcom_split( et2[1] + pjd[1], &pjd[0] );
    et2[0] += pjd[0];

    if ( et2[0] + et2[1] < header->ss[0] || et2[0] + et2[1] > header->ss[1] )
    {
        // fprintf(stderr,"Time is outside ephemeris range.\n");
        retval = -1;
    }
    else
    {
        // Days from start of file.  First part of this calculation should be
        // exact because both values are exactly half integral.
        filetime = ( et2[0] - header->ss[0] ) + et2[1];
        // Data block in file (offset by two in ephcom_readbinary_block to skip two
        // header blocks).
        // blocknum is initially calculated using the convention that the time is
        // in the semi-open interval [datablock[0], datablock[1]), i.e., is
        // strictly less than datablock[1].
        blocknum = (int) ( filetime / header->ss[2] );
        // If time corresponds to datablock[1] of the last block =
        // header->ss[1], then change the above convention to allow
        // time to correspond to that value, that is calculate
        // blocknum index as if the time was just slightly less than
        // header->ss[1].
        if ( et2[0] == header->ss[1] && et2[1] == 0 )
            blocknum--;

        // Read the data block that contains the coefficients for the
        // desired date.
        if ( ephcom_readbinary_block( infp, header, blocknum, datablock ) <= 0 )
        {
            retval = -1;
        }
        else
        {
            // Now step through the bodies and interpolate positions
            // and velocities.

            // Days from block start.  blocktime[0] calculation should
            // be exact because both values exactly half integral.
            blocktime[0] = et2[0] - datablock[0];
            blocktime[1] = et2[1];
            for ( i = 0; i < 13; i++ )
            {
                // The i index corresponds to positions and velocities
                // of solar system bodies, two nutation angles and
                // their time derivatives, or 3 lunar libration angles
                // and their time derivatives as noted in the
                // following table.  (The positions and velocities are
                // solar system barycentric unless noted otherwise.)
                // 0 = Mercury
                // 1 = Venus
                // 2 = Earth-Moon barycenter
                // 3 = Mars
                // 4 = Jupiter
                // 5 = Saturn
                // 6 = Uranus
                // 7 = Neptune
                // 8 = Pluto
                // 9 = Moon (Geocentric)
                // 10 = Sun
                // 11 = nutation angles
                // 12 = lunar librations
                // subspan is always an integer; header->ss[2] is
                // either 2^5 or 2^6 while header->ipt[i][2] is always
                // a low (< 5) power of two.
                subspan     = header->ss[2] / header->ipt[i][2]; // Days/subinterval
                norm_span   = coords->seconds ? subspan * 86400.0 : subspan;
                subinterval = (int) ( ( ( et2[0] - datablock[0] ) + et2[1] ) / subspan );
                // For this corner case calculate the subinterval value as
                // if the time were slightly less than header->ss[1] (which
                // is equal to datablock[1] in this special case).
                if ( et2[0] == header->ss[1] && et2[1] == 0 )
                    subinterval--;

                ncoords    = i == 11 ? 2 : 3; // 2 coords for nutation, else 3
                dataoffset = header->ipt[i][0] - 1 +
                             ncoords * header->ipt[i][1] * subinterval;

                // header->ss[2] / header->ipt[i][2] is always an
                // integer; header->ss[2] is either 2^5 or 2^6 while
                // header->ipt[i][2] is always a low (< 5) power of
                // two.
                subtime = ( blocktime[0] - subinterval * header->ss[2] / header->ipt[i][2] ) + blocktime[1];
                //
                // Divide days in this subblock by total days in
                // subblock to get interval [0,1].  The right part of
                // the expression will evaluate to a whole number:
                // subinterval lengths are all integer multiples of
                // days in a block (all powers of 2).
                //
                chebytime = subtime / subspan;
                chebytime = ( chebytime + chebytime ) - 1.0;
                if ( chebytime < -1.0 || chebytime > 1.0 )
                {
                    fprintf( stderr, "Chebyshev time is beyond [-1,1] interval.\n" );
                    fprintf( stderr,
                        "filetime=%f, blocktime[0]=%f, blocktime[1]=%f, subtime=%f, chebytime=%f\n",
                        filetime, blocktime[0], blocktime[1], subtime, chebytime );
                }
                else
                {
                    //
                    // Everything is as expected.  Interpolate
                    // coefficients to calculate position and velocity
                    // (or angles and angles' time derivatives for the
                    // case of nutation or libration) for the ith
                    // "body" in the solar system at time equivalent
                    // to chebytime.  The number of coordinates is
                    // ncoords which is 3 for everything but nutation
                    // where it is 2.
                    //
                    ephcom_cheby( header->maxcheby, chebytime, norm_span,
                        &datablock[dataoffset],
                        ncoords, header->ipt[i][1], coords->pv[i] );
                }
            }
            //
            //    Set any user-defined coordinates to zero.
            //
            // for (i = 16; i < EPHCOM_NUMOBJECTS; i++)
            //    coords->pv[i][0] = coords->pv[i][1] = coords->pv[i][1] =
            //       coords->pv[i][1] = coords->pv[i][1] = coords->pv[i][1] =  0.0;
            //
            // With interpolations complete, calculate Earth from EMBary and
            // Sun from SSBary.  Preserve other coordinates.
            // N.B. last two elements of nutation are undefined, but
            // as a result of this next loop those locations are zeroed
            // so that all 6 components of each coords->pv[] vector
            // are initialized.
            for ( j = 0; j < 6; j++ )
            {
                // Save original lunar geocentric coords
                coords->pv[15][j] = coords->pv[ 9][j];
                // Save Librations if on file
                coords->pv[14][j] = coords->pv[12][j];
                // Save Nutations if on file.  Last two components
                // are uninitialized so avoid them.
                if ( j < 4 )
                {
                    coords->pv[13][j] = coords->pv[11][j];
                }
                // Save Earth-Moon barycenter.
                coords->pv[12][j] = coords->pv[2][j];
                // Prepare new solar system barycenter coordinates (relative
                // to that center).  Note, this action initializes the
                // last two components of pv[11] (previously used
                // for nutation) for the first time.
                coords->pv[11][j] = 0.;
                //
                // Calculate Earth and Moon from EMBary and geocentric Moon.
                //
                // New Earth
                coords->pv[2][j] -= coords->pv[9][j] / ( 1.0 + header->emrat );
                // New Moon
                coords->pv[9][j] += coords->pv[2][j];
                // The first index corresponds to positions and velocities
                // of solar system bodies, two nutation angles and
                // their time derivatives, or 3 lunar libration angles
                // and their time derivatives as noted in the
                // following table.  (The positions and velocities are
                // solar system barycentric unless noted otherwise. An
                // asterisk preceding a number indicates a change or
                // new index compared to the previous i indices.)
                //  0 = Mercury
                //  1 = Venus
                // *2 = Earth
                //  3 = Mars
                //  4 = Jupiter
                //  5 = Saturn
                //  6 = Uranus
                //  7 = Neptune
                //  8 = Pluto
                // *9 = Moon
                //  10 = Sun
                // *11 = Solar-System barycenter
                // *12 = Earth-Moon barycenter
                // *13 = nutation angles
                // *14 = lunar librations
                // *15 = Moon (geocentric)
            }
            //
            // If we want something other than coordinates relative to
            // the solar system barycenter, subtract coordinates of
            // the reference body (supplied by the calling routine via
            // the coords->center index in the new indexing scheme)
            // from all coordinates except nutation angles (which as a
            // side effect avoids dealing with the 4 components in
            // that special case), libration angles, and geocentric
            // lunar position
            //
            if ( !coords->bary )
            {
                if ( 0 <= coords->center && coords->center <= 12 )
                {
                    for ( i = 0; i < 13; i++ )
                    {
                        if ( i != coords->center )
                        {
                            for ( j = 0; j < 6; j++ )
                                coords->pv[i][j] -= coords->pv[coords->center][j];
                        }
                        else
                        {
                            for ( j = 0; j < 6; j++ )
                                coords->pv[coords->center][j] = 0.;
                        }
                    }
                }
                else
                {
                    fprintf( stderr, "ephcom_get_coords: coords->center = %d is outside the valid range from 0 to 12.\n", coords->center );
                    exit( EXIT_FAILURE );
                }
            }
            if ( !coords->km ) // Calculate AU, not kilometers
            {
                for ( i = 0; i < 15; i++ )
                {
                    // Skip over nutations (which as a side effect
                    // avoids dealing with the 4 components in that
                    // special case) and librations.
                    if ( i == 13 )
                        i = 15;
                    for ( j = 0; j < 6; j++ )
                        coords->pv[i][j] /= header->au;
                }
            }
        }
    }

    return ( retval );
}

//! Interpolate position and velocity (or nutation or libration angles
//! and their time derivatives) at a time point (converted to
//! Chebyshev coordinate in range [-1,1]) using JPL Chebyshev
//! coefficients supplied for one solar system JPL "body" index, where
//! the JPL "bodies" are identified as follows:
//! <table border>
//!   <tr> <td><b>Index</b></td> <td><b>Identification</b></td> </tr>
//!   <tr> <td>0</td> <td>Mercury</td> </tr>
//!   <tr> <td>1</td> <td>Venus</td> </tr>
//!   <tr> <td>2</td> <td>Earth-Moon Barycenter</td> </tr>
//!   <tr> <td>3</td> <td>Mars</td> </tr>
//!   <tr> <td>4</td> <td>Jupiter</td> </tr>
//!   <tr> <td>5</td> <td>Saturn</td> </tr>
//!   <tr> <td>6</td> <td>Uranus</td> </tr>
//!   <tr> <td>7</td> <td>Neptune</td> </tr>
//!   <tr> <td>8</td> <td>Pluto</td> </tr>
//!   <tr> <td>9</td> <td>Moon (Geocentric)</td> </tr>
//!   <tr> <td>10</td> <td>Sun</td> </tr>
//!   <tr> <td>11</td> <td>Nutation angles</td> </tr>
//!   <tr> <td>12</td> <td>Lunar Libration angles</td> </tr>
//! </table>
//!
//! @param maxcoeffs [INPUT ONLY]Maximum number of Chebyshev
//! components possible.
//! @param x [INPUT ONLY]Value of x over [-1,1] for Chebyshev
//! interpolation.
//! @param span [INPUT ONLY]Span of subinterval in the time coordinate
//! used for the time derivatives (velocity or radians per second for
//! the angular coordinates).
//! @param y [INPUT ONLY]Pointer to an array of required Chebyshev
//! coefficients for a particular JPL "body".
//! @param ncoords [INPUT ONLY]Total number of coordinates to
//! interpolate for a particular JPL "body".  This quantity is 3
//! except for nutation where it is two.
//! @param ncoeffs [INPUT ONLY]Number of Chebyshev coefficients per
//! coordinate.
//! @param pv [OUTPUT ONLY]Pointer to an array to hold interpolated
//! positions (or angles) in 1st half, interpolated velocity (or angle
//! time derivatives) in 2nd half for a particular JPL "body".
//!
static inline void
ephcom_cheby( int maxcoeffs, double x, double span, double *y,
              int ncoords, int ncoeffs, double *pv )
{
    int           i, j;
    static double twox;
    static double *pc, *vc;    // Position and velocity polynomial coefficients.
    static double lastx = 2.0; // x from last call; initialize to impossible value
    static int    init  = 1;   // Need to initialize pc[] and vc[]

    //
    //   Allocate position and velocity Chebyshev coefficients.
    //
    if ( init )
    {
        // It is extremely convenient to "permanently" malloc space for
        // pc and vc like this.  Ideally, one would have a special call
        // to ephcom_cheby to free pc and vc or else malloc and free the
        // space outside the routine, but these changes are more trouble
        // than they are worth so we will have to pay the price of
        // valgrind complaining about that unfreed space.
        if ( ( pc = (double *) malloc( maxcoeffs * sizeof ( double ) ) ) == NULL )
        {
            fprintf( stderr, "ephcom_cheby: Cannot malloc pc" );
            exit( EXIT_FAILURE );
        }
        if ( ( vc = (double *) malloc( maxcoeffs * sizeof ( double ) ) ) == NULL )
        {
            fprintf( stderr, "ephcom_cheby: Cannot malloc vc" );
            exit( EXIT_FAILURE );
        }
        init = 0;
    }
    //
    //   This need only be called once for each Julian Date,
    //   saving a lot of time initializing polynomial coefficients.
    //
    if ( lastx != x )
    {
        lastx = x;
        twox  = x + x; // For Chebyshev recursion
        //
        //           Initialize position polynomial coefficients
        //
        pc[0] = 1.0;   // Chebyshev T[0](x) = 1
        pc[1] = x;     // Chebyshev T[1](x) = x
        for ( i = 2; i < maxcoeffs; i++ )
        {
            pc[i] = twox * pc[i - 1] - pc[i - 2];
            // Resolve bug with gcc generating -0.0 (also makes
            // the smallest represented number equal to zero).
            //
            if ( pc[i] * pc[i] == 0.0 )
            {
                pc[i] = 0.0;
            }
        }
        //
        // Initialize derivative polynomial coefficients
        //
        vc[0] = 0.0;         // d(1)/dx        = 0
        vc[1] = 1.0;         // d(x)/dx        = 1
        vc[2] = twox + twox; // d(2x^2 - 1)/dx = 4x
        for ( i = 3; i < maxcoeffs; i++ )
        {
            vc[i] = twox * vc[i - 1] + pc[i - 1] + pc[i - 1] - vc[i - 2];
        }
    }
    //
    //   Interpolate to get position for each component
    //
    for ( i = 0; i < ncoords; i++ ) // Once each for x, y, and z
    {
        pv[i] = 0.0;
        for ( j = ncoeffs - 1; j >= 0; j-- )
        {
            pv[i] += pc[j] * y[i * ncoeffs + j];
        }
    }
    //
    // Interpolate velocity (first derivative)
    //
    for ( i = 0; i < ncoords; i++ )
    {
        pv[ncoords + i] = 0.0;
        for ( j = ncoeffs - 1; j >= 0; j-- )
        {
            pv[ncoords + i] += vc[j] * y[i * ncoeffs + j];
        }
        pv[ncoords + i] *= 2.0 / span;
    }
}

//! Convert Julian Day to calendar date and time.  From pp. 604, 606
//! in the Explanatory Supplement to the Astronomical Almanac.
//!
//! @param tjd [IN ONLY]Double-precision Julian Day number.
//! @param idate [OUT ONLY] integer array of 6 values which upon
//! return will contain the integer year, month, day, hour, minute,
//! and second corresponding to tjd.
//! @param calendar_type [IN ONLY]Integer value which controls the
//! kind of calendar used for the transformation: -1=Julian;
//! 0=automatic; 1=Gregorian.  If automatic, use Julian calendar for
//! dates before 15 October 1582.
//!
static void
ephcom_jd2cal( double tjd, int idate[6], int calendar_type )
{
    int ihour, imin, isec;
    int j;
// From Explanatory Supplement to Astronomical Almanac, pp. 604, 606
    int I, J, K, L, N, D, M, Y;

    tjd  += 0.5 + 0.5 / 86400.0; // Round to nearest second
    j     = tjd;                 // Integer Julian Day
    tjd   = ( tjd - j ) * 24.0;
    ihour = tjd;
    tjd   = ( tjd - ihour ) * 60.0;
    imin  = tjd;
    tjd   = ( tjd - imin ) * 60.0;
    isec  = tjd;
//
//   Julian calendar.  Explanatory Supplement to Astronomical Alamanac, p. 606.
//   If automatic, use Julian calendar for dates before 15 October 1582.
//
    if ( calendar_type == -1 || ( calendar_type == 0 && j <= 2299160 ) )
    {
        J = j + 1402;
        K = ( J - 1 ) / 1461;
        L = J - 1461 * K;
        N = ( L - 1 ) / 365 - L / 1461;
        I = L - 365 * N + 30;
        J = ( 80 * I ) / 2447;
        D = I - ( 2447 * J ) / 80;
        I = J / 11;
        M = J + 2 - 12 * I;
        Y = 4 * K + N + I - 4716;
    }