ephcom.c

//! @file
//! Source code for the ephcom library.
//!

// Copyright (C) 1994-2004 Paul Hardy
// Copyright (C) 2011 Alan W. Irwin
//
// This file is part of the timeephem software project.
//
// timeephem is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published
// by the Free Software Foundation; version 2 of the License.
//
// timeephem is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with timeephem; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA

//   Note that this software is not a product of the Jet Propulsion
//   Laboratory; it just uses and supports their ASCII and binary
//   ephemeris files.  Please don't mail JPL concerning any bugs.
//   Send bug reports or suggestions to airwin@users.sourceforge.net instead.
//
//   This file contains the following routines.  Running "make test" will
//   check the proper operation of every one of these routines.
//
//      ephcom_readascii_header()   - read ASCII header file
//      ephcom_readascii_block()    - read ASCII coefficient block
//      ephcom_readbinary_header()  - read header from binary ephemeris
//      ephcom_readbinary_block()   - read coefficient block from binary
//                                    ephemeris
//      ephcom_writeascii_header()  - write header in ASCII
//      ephcom_writeascii_block()   - write coefficient block in ASCII
//      ephcom_writebinary_header() - write header to binary ephemeris file
//      ephcom_writebinary_block()  - write coefficient block to binary
//                                    ephemeris file
//      ephcom_parse_block()        - parse ("pretty print") a coefficient block
//      ephcom_nxtgrp()             - read next "GROUP" from ASCII header
//      ephcom_outdouble()          - write byte-swapped double to a file
//      ephcom_outint()             - write byte-swapped int to a file
//      ephcom_indouble()           - read byte-swapped double from a file
//      ephcom_inint()              - read byte-swapped int from a file
//      ephcom_doublstrc2f()        - change C ASCII double string to FORTRAN
//                                    [there is no corresponding
//                                     ephcom_doublestrf2c() routine;
//                                     for FORTRAN to C format conversion,
//                                     just change FORTRAN's double precision
//                                     'D' exponent to 'E' in your software
//                                     and everything else should parse fine]
//      ephcom_pleph()              - calculate <x,y,z> and <xdot,ydot,zdot>
//                                    for a given target and center, AFTER
//                                    calling ephcom_get_coords() (different
//                                    sequence than with JPL's FORTRAN PLEPH)
//      ephcom_get_coords()         - calculate <x,y,z> and <xdot,ydot,zdot>
//                                    for all Solar System objects at a
//                                    given time
//      ephcom_cheby()              - interpolates Chebyshev coefficients
//                                    for one sub-block of coefficients
//                                    for one Solar System object at a
//                                    given time
//      ephcom_jd2cal()             - convert Julian Day to Julian or Gregorian
//                                    Year, Month, Day, Hour, Minute, Second
//      ephcom_cal2jd()             - convert Julian or Gregorian calendar
//                                    Year, Month, Day, Hour, Minute, Second
//                                    to Julian Day
//
//   The indouble() and outdouble() routines rely upon the gnulliver64c()
//   routine from gnulliver.c.
//
//   The inint() and outint() routines rely upon the gnulliver32() routine
//   from gnulliver.c.
//
//
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <math.h> // IF_SAME_DATE macro uses fabs, ephcom_split uses modf.
#include "ephcom.h"
#include "gnulliver.h"

// Declare static functions.
static void
ephcom_nxtgrp( char *group, const char *expected, FILE *infile );

static void
ephcom_outdouble( FILE *outfp, double x );

static void
ephcom_outint( FILE * outfp, unsigned u );

static void
ephcom_doublestrc2f( char *buf );

static double
ephcom_exact_time( double time );

static double
ephcom_split( double time, double * itime );

// ephcom_cheby() - interpolate at a point using Chebyshev coefficients.
static inline void
ephcom_cheby( int maxcoeffs, double x, double span, double *y,
              int ncoords, int ncoeffs, double *pv );

//
//   ephcom_jd2cal() - convert Julian Day to calendar date and time.
//
//      tjd: double precision Julian Day
//      idate: integer year, month, day, hour, minute, second of tjd
//      calendar_type: -1=Julian; 0=Automatic; 1=Gregorian
//
//   If automatic, use Julian calendar for dates before 15 October 1582.
//
//   From pp. 604, 606 in the Explanatory Supplement to the Astronomical Almanac.
//
static void
ephcom_jd2cal( double tjd, int idate[6], int calendar_type );

// Start of function definitions.

//! Read a JPL ephemeris ASCII header from the FILE pointed to by the
//! infp argument and store all header data in the ephcom_Header
//! struct pointed to by the header argument.  If any errors are
//! detected this routine writes a message to stderr and exits.
//!
//! @param infp [IN ONLY]Pointer to an ascii version of a JPL
//! ephemeris FILE.
//! @param header [OUT ONLY]Pointer to an ephcom_Header struct which
//! upon return will contain all the JPL ephemeris header information
//! that has been read from the FILE.
//!
void
ephcom_readascii_header( FILE * infp, ephcom_Header *header )
{
    char   group[13];        // To store the "GROUP" header line information
    double val1, val2, val3; // To read text line with 3 double precision words
    int    i, j, k, n;
    int    iword;            // word number we're reading in a line
    int    blockout;         // number of bytes we've written to current block/rec in file
    int    blockbytes;       // number of bytes in a block, equals 8 * ncoeff

    char   readbuf[EPHCOM_MAXLINE + 1];

    char   outhcars[EPHCOM_MAXLINE + 1];
    size_t fwrite( const void *ptr, size_t size, size_t nmemb, FILE *stream );

//
//   First header line: KSIZE= # NCOEFF= #
//
    if ( infp != stdin )
        rewind( infp );
    fgets( readbuf, EPHCOM_MAXLINE, infp );
    sscanf( readbuf, "%*6s%6d%*11s%6d", &header->ksize, &header->ncoeff );
    blockbytes = 8 * header->ncoeff; // The size of a double, times # of doubles/block
    if ( header->ksize != 2 * header->ncoeff )
    {
        fprintf( stderr, "Badly formed header; header->ksize != 2*header->ncoeff\n\n" );
        exit( 1 );
    }
//
//   GROUP 1010: Title of ephemeris (DE/LE number, start JD, end JD)
//
//
//   Blank all of header->ttl.  Note that three fgets below
//   only defines part of ttl so this blanking keeps valgrind
//   quiet for subsequent accesses to all of ttl.
//
    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < EPHCOM_MAXTTL; j++ )
            header->ttl[i][j] = ' ';
        header->ttl[i][EPHCOM_MAXTTL] = '\0';
    }
    ephcom_nxtgrp( group, "GROUP   1010", infp );
    fgets( header->ttl[0], EPHCOM_MAXTTL + 2, infp ); // JPL Ephemeris title line
    if ( strncmp( header->ttl[0], "JPL ", 4 ) != 0 )
    {
        fprintf( stderr, "\nERROR: file is not a JPL ASCII header file.\n\n" );
        exit( 1 );
    }
    fgets( header->ttl[1], EPHCOM_MAXTTL + 2, infp ); // Start epoch
    fgets( header->ttl[2], EPHCOM_MAXTTL + 2, infp ); // Finish epoch
//
//   Convert any newlines or tabs to single spaces.
//
    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < EPHCOM_MAXTTL; j++ )
            if ( isspace( header->ttl[i][j] ) )
                header->ttl[i][j] = ' ';
        header->ttl[i][EPHCOM_MAXTTL] = '\0';
    }
//
//   GROUP 1030: Start and End JD, timestep (in JD) per block.
//
    ephcom_nxtgrp( group, "GROUP   1030", infp );
    fgets( readbuf, EPHCOM_MAXLINE, infp );
    sscanf( readbuf, " %lE %lE %lE", &header->ss[0], &header->ss[1], &header->ss[2] );
//
//   GROUP 1040: Constant names.
//
    ephcom_nxtgrp( group, "GROUP   1040", infp );
    fgets( readbuf, EPHCOM_MAXLINE, infp );
    header->ncon = atoi( readbuf );
//
//   Now read the constant names, 10 per line, each 6 characters long
//   preceded by 2 blanks.  Pad names with blanks to make 6 characters.
//
    for ( i = 0; i < header->ncon; )
    {
        fgets( readbuf, EPHCOM_MAXLINE, infp );
        if ( ( j = strlen( readbuf ) ) < 81 ) // Pad end with blanks for copying
        {
            // initial j is such a value that readbuf[j-1] is '\n'
            while ( j < 81 )
                readbuf[j++ - 1] = ' ';
            readbuf[80] = '\n';
            readbuf[81] = '\0';
        }
        for ( iword = 0; iword < 10 && i < header->ncon; iword++, i++ )
        {
            strncpy( header->cnam[i], &readbuf[2 + iword * 8], 6 );
            header->cnam[i][6] = '\0';
        }
    }
//
//   GROUP 1041: Constant values.
//
    ephcom_nxtgrp( group, "GROUP   1041", infp );
    fgets( readbuf, EPHCOM_MAXLINE, infp );
    header->nval = atoi( readbuf );
    if ( header->nval != header->ncon )
    {
        fprintf( stderr, "Error: number of constants and values not equal.\n\n" );
        exit( 1 );
    }
//
//   Now read constant values, 3 per line, 26 characters each.
//
    for ( i = 0; i < header->ncon; i += 3 )
    {
        fgets( readbuf, EPHCOM_MAXLINE, infp );
        for ( j = 0; j < strlen( readbuf ); j++ )
            if ( tolower( readbuf[j] ) == 'd' )
                readbuf[j] = 'E';
        // exponent is 'E'
        sscanf( readbuf, "%lE %lE %lE",
            &header->cval[i], &header->cval[i + 1], &header->cval[i + 2] );
    }
//
//   GROUP 1050: Constant values.
//
    ephcom_nxtgrp( group, "GROUP   1050", infp );
    for ( i = 0; i < 3; i++ )
    {
        fgets( readbuf, EPHCOM_MAXLINE, infp ); // Read line of 13 6-digit integers
        for ( j = 0; j < 12; j++ )
        {
            header->ipt[j][i] = atoi( &readbuf[6 * j] );
        }
        header->lpt[i] = atoi( &readbuf[6 * 12] );
    }
//
//   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;
//
//   Set the maximum number of Chebyshev coefficients possible for this file,
//   to initialize position and velocity Chebyshev coefficient arrays during
//   Chebyshev interpolation.
//
    header->maxcheby = 0;
    for ( i = 0; i < 12; i++ )
        if ( header->ipt[i][1] > header->maxcheby )
            header->maxcheby = header->ipt[i][1];
    if ( header->lpt[1] > header->maxcheby )
        header->maxcheby = header->lpt[1];

    header->au    = 0.0;
    header->emrat = 0.0;
    header->numde = 0;
    for ( i = 0; i < header->ncon; i++ )
    {
        if ( strncmp( header->cnam[i], "AU    ", 6 ) == 0 )
            header->au = header->cval[i];
        else if ( strncmp( header->cnam[i], "EMRAT ", 6 ) == 0 )
            header->emrat = header->cval[i];
        else if ( strncmp( header->cnam[i], "DENUM ", 6 ) == 0 )
            header->numde = header->cval[i];
        else if ( strncmp( header->cnam[i], "CLIGHT", 6 ) == 0 )
            header->clight = header->cval[i];
        else if ( strncmp( header->cnam[i], "LENUM ", 6 ) == 0 )
            header->numle = header->cval[i];
    }
    if ( header->numle == 0 )
        header->numle = header->numde;
//
//   GROUP 1070: Constant values.
//
    ephcom_nxtgrp( group, "GROUP   1070", infp );
//
//   Now we're pointing to the first block of coefficient data, after header.
//   Return at the point where we can start reading coefficients.
//
}

//! Read a block of data coefficients from a JPL ASCII ephemeris file.
//!
//! @param infp [IN ONLY]Pointer to an ascii 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 the FILE.
//! @param datablock [OUT ONLY]Pointer to an array that upon
//! successful return will be filled with datapoints == header->ncoeff
//! data points.
//! @returns number of coefficients read or 0 at EOF or some other i/o
//! error.
//!
int
ephcom_readascii_block(
    FILE * infp,
    ephcom_Header *header,
    double *datablock )
{
    int    i, j;
    int    datalines;        // lines of data we've read
    int    datapoints;       // points of data we've read/converted/written
    char   readbuf[EPHCOM_MAXLINE + 1];
    double val1, val2, val3; // To read text line with 3 double precision words

//
//   First line in an ASCII block will be the block number, followed by
//   the number of coefficients.
//
    datalines  = 0; // Not reported, but leave in for debugging
    datapoints = 0;
    if ( fgets( readbuf, EPHCOM_MAXLINE, infp ) && !feof( infp ) )
    {
        sscanf( readbuf, "%d %d", &i, &j );
        if ( j != header->ncoeff )
        {
            fprintf( stderr,
                "\nERROR: ASCII data file's %d coefficients/block\n", j );
            fprintf( stderr,
                "       doesn't match ASCII header's %d coefficients/block.\n\n",
                header->ncoeff );
            exit( 1 );
        }
        datalines++;
        while ( datapoints < header->ncoeff && !feof( infp ) )
        {
            fgets( readbuf, EPHCOM_MAXLINE, infp );
            for ( j = 0; j < strlen( readbuf ); j++ )
                if ( tolower( readbuf[j] ) == 'd' )
                    readbuf[j] = 'e';
            datalines++;
            //
            // This is horrible, but use "%le" here and "%lE in the other
            // ASCII data routine (ephcom_readascii_header) so gcc won't try
            // to store the formats in the same location and write to them.
            // (Problem with gcc not acting like K&R without -traditional flag
            // and without -fwritable-strings flag.)
            //
            sscanf( readbuf, " %le %le %le", &val1, &val2, &val3 );
            datablock[datapoints++] = val1;
            if ( ( datapoints ) < header->ncoeff )
            {
                datablock[datapoints++] = val2;
                if ( datapoints < header->ncoeff )
                {
                    datablock[datapoints++] = val3;
                }
            }
        }
    }
    return ( datapoints );
}

//! Read a JPL ephemeris binary header from the FILE pointed to by the
//! infp argument and store all header data in the ephcom_Header
//! struct pointed to by the header argument.  If any errors are
//! detected this routine writes a message to stderr and exits.
//!
//! @param infp [IN ONLY]Pointer to a binary version of a JPL
//! ephemeris FILE.
//! @param header [OUT ONLY]Pointer to an ephcom_Header struct which
//! upon return will contain all the JPL ephemeris header information
//! that has been read from the FILE.
//!
void
ephcom_readbinary_header( FILE * infp, ephcom_Header *header )
{
    int i, j, k;

    if ( infp != stdin )
        rewind( infp );
//
//   Read title lines.
//
    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < EPHCOM_MAXTTL; j++ )
        {
            header->ttl[i][j] = fgetc( infp );
        }
        if ( i == 0 && strncmp( header->ttl[0], "JPL ", 4 ) != 0 )
        {
            fprintf( stderr, "\nERROR: file is not a JPL ephemeris file.\n\n" );
            if ( strncmp( header->ttl[0], "KSIZE", 5 ) == 0 )
                fprintf( stderr, "File is an ASCII JPL ephemeris header instead.\n\n" );
            exit( 1 );
        }
        header->ttl[i][j] = '\0';
    }
//
//   Read constant names.
//
    for ( i = 0; i < 400; i++ )
    {
        for ( j = 0; j < 6; j++ )
        {
            header->cnam[i][j] = fgetc( infp );
        }
        header->cnam[i][j] = '\0';
    }
//
//   Read ephemeris start epoch, stop epoch, and step size (in Julian Days).
//
    for ( i = 0; i < 3; i++ )
    {
        header->ss[i] = ephcom_indouble( infp );
    }
    // These values are half integral Julian dates.  Make sure there is no
    // numerical noise in these values.
    header->ss[0] = ephcom_exact_time( header->ss[0] );
    header->ss[1] = ephcom_exact_time( header->ss[1] );
    // This value is an integral number of days (a power of two).  Make sure there
    // is no numerical noise in this value.
    header->ss[2] = (double) (int) ( header->ss[2] + 0.01 );
//
//   Read NCON, AU, EMRAT.
//
    header->ncon  = ephcom_inint( infp );
    header->au    = ephcom_indouble( infp );
    header->emrat = ephcom_indouble( infp );
    header->nval  = header->ncon;
//
//   Read indexes for coefficients in data block.  Written in transposed
//   order (Fortran and C matrices are transposed).
//
    for ( i = 0; i < 12; i++ )
    {
        for ( j = 0; j < 3; j++ )
        {
            header->ipt[i][j] = ephcom_inint( infp );
        }
    }
    header->numde = ephcom_inint( infp ); // Get ephemeris number
    for ( i = 0; i < 3; i++ )
        header->lpt[i] = ephcom_inint( infp );
//
//   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;
//
//   Set the maximum number of Chebyshev coefficients possible for this file,
//   to initialize position and velocity Chebyshev coefficient arrays during
//   Chebyshev interpolation.
//
    header->maxcheby = 0;
    for ( i = 0; i < 12; i++ )
        if ( header->ipt[i][1] > header->maxcheby )
            header->maxcheby = header->ipt[i][1];
    if ( header->lpt[1] > header->maxcheby )
        header->maxcheby = header->lpt[1];

//
//    From JPL ephemeris number, set NCOEFF and calculate KSIZE = 2*NCOEFF.
//
// switch (header->numde) {
//    case 102:
//       header->ncoeff = 773;
//       break;
//    case 200:
//       header->ncoeff = 826;
//       break;
//    case 202:
//       header->ncoeff = 826;
//       break;
//    case 403:
//       header->ncoeff = 1018;
//       break;
//    case 405:
//       header->ncoeff = 1018;
//       break;
//    case 406:
//       header->ncoeff = 728;
//       break;
//    case 410:
//       header->ncoeff = 1018;
//       break;
//    default:
//       header->ncoeff = 1018;
//       break;
//    }
//
//   Calculate number of coefficients, starting with
//   highest index into a data block (stored in j).
//
    j = 0;
    for ( i = 1; i < 12; i++ )
        if ( header->ipt[i][1] > 0 && header->ipt[i][0] > header->ipt[j][0] )
            j = i;

//
//   Now see if the starting point we found is lower than where
//   lpt[] starts.  If not, use lpt[] for largest value.
//
    if ( header->lpt[1] > 0 && header->lpt[0] > header->ipt[j][0] )
    {
        header->ncoeff = header->lpt[0] - 1 + // starting point
                         ( header->lpt[1] *   // coefficients per coordinate
                           header->lpt[2] ) * // subblocks per block
                         3;                   // coordinates
    }
    else
    {
        header->ncoeff = header->ipt[j][0] - 1 + // starting point
                         ( header->ipt[j][1] *   // coefficients per coordinate
                           header->ipt[j][2] ) * // subblocks per block
                         ( j == 11 ? 2 : 3 );    // coordinates
    }

    header->ksize = header->ncoeff + header->ncoeff; // KSIZE = 2*NCOEFF always
//
//   Skip to second block in file.
//
    fseek( infp, header->ncoeff * 8, SEEK_SET );
//
//   Read ephemeris constants.
//
    for ( i = 0; i < header->ncon; i++ )
    {
        header->cval[i] = ephcom_indouble( infp );
        if ( strncmp( header->cnam[i], "NCOEFF", 6 ) == 0 )
        {
            header->ncoeff = header->cval[i];
            header->ksize  = 2 * header->cval[i];
        }
        else if ( strncmp( header->cnam[i], "LENUM ", 6 ) == 0 )
            header->numle = header->cval[i];
    }
    if ( header->numle == 0 )
        header->numle = header->numde;
}

//! Read a block of data coefficients from a JPL binary ephemeris file.
//!
//! @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 the infp FILE.
//! @param blocknum [IN ONLY]Requested direct access block number with
//! zero corresponding to the first data block after the two header
//! blocks in the binary ephemeris.
//! @param datablock [OUT ONLY]Pointer to an array that upon
//! successful return will be filled with datapoints == header->ncoeff
//! data points.
//! @returns number of coefficients read or 0 at EOF or some other i/o
//! error.
//!
int
ephcom_readbinary_block(
    FILE *infp,
    ephcom_Header *header,
    int blocknum,
    double *datablock
    )
{
    int  i;
    long filebyte;

    filebyte = ( blocknum + 2 ) * header->ncoeff * 8; // 8 bytes per coefficient
    fseek( infp, filebyte, SEEK_SET );
    for ( i = 0; !feof( infp ) && !ferror( infp ) && i < header->ncoeff; i++ )
    {
        datablock[i] = ephcom_indouble( infp );
    }

    if ( feof( infp ) || ferror( infp ) )
        i = 0;    // 0 --> EOF or any other i/o error.

    // First two values of data block are half-integral Julian dates.  Make sure
    // there is no numerical noise in these data.  (This makes a difference for
    // early versions of de422.)
    if ( i >= 1 )
    {
        datablock[0] = ephcom_exact_time( datablock[0] );
        datablock[1] = ephcom_exact_time( datablock[1] );
    }

    return ( i ); // Number of coefficients successfuly read (all or nothing).
}

//! Write JPL ephemeris header information in ASCII format.
//! @param outfp [IN ONLY]Pointer to the FILE which will be used to
//! output the formatted header information.
//! @param header [IN and OUT]Pointer to an ephcom_Header struct which
//! contains the JPL ephemeris header information to be output to the
//! FILE.  However, some self-consistency adjustments of the header
//! are made before such FILE output and those are returned to the
//! calling routine as well.
//!
void
ephcom_writeascii_header( FILE * outfp, ephcom_Header *header )
{
    char        group[13];
    double      val1, val2, val3; // To read text line with 3 double precision words
    int         i, j, k, n;
    int         iword;            // word number we're reading in a line
    int         blockout;         // number of bytes we've written to current block/rec in file
    int         blockbytes;       // number of bytes in a block, equals 8 * ncoeff
    static char spaces[EPHCOM_MAXTTL] = "                                                                                 \n";
    int         idate[6];
    const char        *month[12] = { "JAN", "FEB", "MAR", "APR", "MAY", "JUN",
                               "JUL", "AUG", "SEP", "OCT", "NOV", "DEC" };

    char        writebuf[EPHCOM_MAXLINE + 1];
    char        outhcars[EPHCOM_MAXLINE + 1];
    size_t fwrite( const void *ptr, size_t size, size_t nmemb, FILE *stream );

//
//   First header line: KSIZE= # NCOEFF= #
//
    blockbytes = 8 * header->ncoeff; // sizeof(double) * # of doubles/block
    fprintf( outfp, "KSIZE=%5d    NCOEFF=%5d\n", header->ksize, header->ncoeff );
    if ( header->ksize != 2 * header->ncoeff )
    {
        fprintf( stderr, "Badly formed header; KSIZE <> 2*NCOEFF\n" );
        exit( 1 );
    }
//
//   GROUP 1010: Title of ephemeris (DE/LE number, start JD, end JD)
//
    fprintf( outfp, " \n" ); // blank line
    fprintf( outfp, "GROUP   1010\n" );
    fprintf( outfp, " \n" ); // blank line
//
//   Header title lines with dates, for example:
//
//      JPL Planetary Ephemeris DE405/LE405
//      Start Epoch: JED=  2305424.5 1599 DEC 09 00:00:00
//      Final Epoch: JED=  2525008.5 2201 FEB 20 00:00:00
//
    sprintf( header->ttl[0], "JPL Planetary Ephemeris DE%03d/LE%03d",
        header->numde, header->numle );
    k = strlen( header->ttl[0] );
    strcpy( &header->ttl[0][k], &spaces[k] );
    ephcom_jd2cal( header->ss[0], idate, 0 );
    sprintf( header->ttl[1], "Start Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
        header->ss[0], idate[0], month[idate[1] - 1], idate[2],
        idate[3], idate[4], idate[5] );
    k = strlen( header->ttl[1] );
    strcpy( &header->ttl[1][k], &spaces[k] );
    ephcom_jd2cal( header->ss[1], idate, 0 );
    sprintf( header->ttl[2], "Final Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
        header->ss[1], idate[0], month[idate[1] - 1], idate[2],
        idate[3], idate[4], idate[5] );
    k = strlen( header->ttl[2] );
    strcpy( &header->ttl[2][k], &spaces[k] );

//
//   Don't print trailing blanks at the end of these 3 lines.
//
    for ( i = 0; i < 3; i++ )
    {
        strncpy( writebuf, header->ttl[i], EPHCOM_MAXTTL + 1 );
        for ( j = EPHCOM_MAXTTL; isspace( writebuf[j] ) || writebuf[j] == '\0'; j-- )
            writebuf[j] = '\0';
        // To match end space in JPL Epoch header lines.
        if ( i > 0 )
            writebuf[++j] = ' ';

        fprintf( outfp, "%s\n", writebuf );
    }
//
//   GROUP 1030: Start and End JD, timestep (in JD) per block.
//
    fprintf( outfp, " \n" ); // blank line
    fprintf( outfp, "GROUP   1030\n" );
    fprintf( outfp, " \n" ); // blank line

    fprintf( outfp, "%12.2f%12.2f%11.0f.\n",
        header->ss[0], header->ss[1], header->ss[2] );
//
//   GROUP 1040: Constant names.
//
    fprintf( outfp, " \n" ); // blank line
    fprintf( outfp, "GROUP   1040\n" );
    fprintf( outfp, " \n" ); // blank line

    fprintf( outfp, "%6d\n", header->ncon );

//
//   Now write the constant names, 10 per line, each 6 characters long
//   preceded by 2 blanks.  Pad names with blanks to make 6 characters.
//
    for ( i = 0; i < header->ncon; i++ )
    {
        fprintf( outfp, "  %-6s", header->cnam[i] );
        if ( i % 10 == 9 )
            fprintf( outfp, "\n" );
    }
    if ( i % 10 != 0 ) // Pad last line with spaces (i is 1 more than above)
    {
        for (; i % 10 != 0; i++ )
            fprintf( outfp, "        " );
        fprintf( outfp, "\n" );
    }
//
//   GROUP 1041: Constant values.
//
    fprintf( outfp, " \n" ); // blank line
    fprintf( outfp, "GROUP   1041\n" );
    fprintf( outfp, " \n" ); // blank line

    fprintf( outfp, "%6d\n", header->nval );

    if ( header->nval != header->ncon )
    {
        fprintf( stderr, "Error: number of constants and values not equal.\n\n" );
        exit( 1 );
    }
//
//   Now write constant values, 3 per line, 26 characters each.
//
    for ( i = 0; i < header->ncon; i += 3 )
    {
        val1 = header->cval[i];
        val2 = ( i + 1 < header->ncon ) ? header->cval[i + 1] : 0.0;
        val3 = ( i + 2 < header->ncon ) ? header->cval[i + 2] : 0.0;

        // Write values, 3 coefficients per line, pad lines with 0.0E+00
        // Must have trailing blank to make room for reformatted
        // fortran version (with leading "0.") created below.
        // Note there is (just) room for Windows 3-digit exponent in
        // the format.
        sprintf( writebuf, "%25.17E %25.17E %25.17E ", val1, val2, val3 );

        // Now re-format numbers the way the JPL header file writes them:
        // all with a leading "0.", so the exponent is one greater.
        // If the number is written in Windows 3-digit exponent format,
        // then it is shifted by one byte to the right to overwrite
        // the assumed leading 0 of the exponent (or error out if
        // there are three exponent digits but the leading one is not
        // zero, i.e., this logic will not work for numbers greater
        // than or equal to 1.e100 or less than or equal to 1.e-101, but
        // this limitation is also in the current JPL ascii format which
        // we are trying to mimic as closely as possible here.
        ephcom_doublestrc2f( &writebuf[0] );  // Reformat first number
        ephcom_doublestrc2f( &writebuf[26] ); // Reformat second number
        ephcom_doublestrc2f( &writebuf[52] ); // Reformat third number
        fprintf( outfp, "%s\n", writebuf );
    }
//
//   GROUP 1050: Constant values.
//
    fprintf( outfp, " \n" ); // blank line
    fprintf( outfp, "GROUP   1050\n" );
    fprintf( outfp, " \n" ); // blank line
//
//   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;
//
//   Write ipt array in transposed order (arrays are transposed in FORTRAN
//   from their order in C).
//
    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < 12; j++ )
        {
            fprintf( outfp, "%6d", header->ipt[j][i] );
        }
        fprintf( outfp, "%6d\n", header->lpt[i] );
    }
//
//   GROUP 1070: Constant values.
//
    fprintf( outfp, " \n" ); // blank line
    fprintf( outfp, "GROUP   1070\n" );
    fprintf( outfp, " \n" ); // blank line
//
//   Now we're pointing to the first block of coefficient data, after header.
//
}

//! Write JPL ephemeris coefficient block of data in ASCII format.
//! @param outfp [IN ONLY]Pointer to the FILE which will be used to
//! output the formatted 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 included in the formatted block and some of these
//! data are used to control how much is written in that formatted block.
//! @param blocknum [IN ONLY]Block number.  This value + 1 will be used
//! in the output formatted block of data.
//! @param datablock [IN ONLY]Pointer to an array that contains the
//! block of data that will be output in formatted form.
//!
void
ephcom_writeascii_block(
    FILE * outfp,
    ephcom_Header *header,
    int blocknum,
    double *datablock )
{
    double val1, val2, val3; // To write text line with 3 double precision words
    int    i, j, k, n;

    char   writebuf[EPHCOM_MAXLINE + 1];
    char   outhcars[EPHCOM_MAXLINE + 1];
    size_t fwrite( const void *ptr, size_t size, size_t nmemb, FILE *stream );
    int fputc( int, FILE * );

//
//   Write first line in block, which is block number and ncoeff.
//
    fprintf( outfp, "%6d%6d\n", blocknum + 1, header->ncoeff );
//
//   Now write the data, 3 coefficients per line, 26 characters each.
//   Convert format to match what appears in JPL Ephemeris ASCII files.
//
    for ( i = 0; i < header->ncoeff; i += 3 )
    {
        val1 = datablock[i];
        val2 = ( i + 1 ) < header->ncoeff ? datablock[i + 1] : 0.0;
        val3 = ( i + 2 ) < header->ncoeff ? datablock[i + 2] : 0.0;

        // Write values, 3 coefficients per line, pad lines with 0.0E+00
        // Must have trailing blank to make room for reformatted
        // fortran version (with leading "0.") created below.
        // Note there is (just) room for Windows 3-digit exponent in
        // the format.
        sprintf( writebuf, "%25.17E %25.17E %25.17E ", val1, val2, val3 );

        // Now re-format numbers the way the JPL header file writes them:
        // all with a leading "0.", so the exponent is one greater.
        // If the number is written in Windows 3-digit exponent format,
        // then it is shifted by one byte to the right to overwrite
        // the assumed leading 0 of the exponent (or error out if
        // there are three exponent digits but the leading one is not
        // zero, i.e., this logic will not work for numbers greater
        // than or equal to 1.e100 or less than or equal to 1.e-101, but
        // this limitation is also in the current JPL ascii format which
        // we are trying to mimic as closely as possible here.
        ephcom_doublestrc2f( &writebuf[0] );  // Reformat first number
        ephcom_doublestrc2f( &writebuf[26] ); // Reformat second number
        ephcom_doublestrc2f( &writebuf[52] ); // Reformat third number
        fprintf( outfp, "%s\n", writebuf );
    }
}

//! Write JPL ephemeris header information in binary form.
//! @param outfp [IN ONLY]Pointer to the FILE which will be used to
//! output the binary header information.
//! @param header [IN ONLY]Pointer to an ephcom_Header struct which
//! contains the JPL ephemeris header information to be output in
//! binary form.
//!
void
ephcom_writebinary_header( FILE *outfp, ephcom_Header *header )
{
    char   readbuf[EPHCOM_MAXLINE + 1];
    char   group[13];        // To hold "GROUP" header line
    int    blockout;         // number of bytes we've written to current block/rec in file
    int    blockbytes;       // number of bytes in a block, equals 8 * ncoeff

    double val1, val2, val3; // To read text line with 3 double precision words
    int    i, j, k, n;
    int    idate[6];
    const char   *month[12] = { "JAN", "FEB", "MAR", "APR", "MAY", "JUN",
                          "JUL", "AUG", "SEP", "OCT", "NOV", "DEC" };

    char   outhcars[EPHCOM_MAXLINE + 1];
    size_t fwrite( const void *ptr, size_t size, size_t nmemb, FILE *stream );

//
//   Point to beginning of output file.
//
    rewind( outfp );
//
//   First header line: KSIZE= # NCOEFF= #
//
    blockbytes = sizeof ( double ) * header->ncoeff;
//
//   Start writing output ephemeris, beginning with header.
//
//
//   Header title lines with dates, for example:
//
//      JPL Planetary Ephemeris DE405/LE405
//      Start Epoch: JED=  2305424.5 1599 DEC 09 00:00:00
//      Final Epoch: JED=  2525008.5 2201 FEB 20 00:00:00
//
    sprintf( header->ttl[0], "JPL Planetary Ephemeris DE%03d/LE%03d",
        header->numde, header->numle );
    for ( i = strlen( header->ttl[0] ); i < EPHCOM_MAXTTL; i++ )
        header->ttl[1][i] = ' ';
    ephcom_jd2cal( header->ss[0], idate, 0 );
    sprintf( header->ttl[1], "Start Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
        header->ss[0], idate[0], month[idate[1] - 1], idate[2],
        idate[3], idate[4], idate[5] );
    for ( i = strlen( header->ttl[1] ); i < EPHCOM_MAXTTL; i++ )
        header->ttl[1][i] = ' ';
    ephcom_jd2cal( header->ss[1], idate, 0 );
    sprintf( header->ttl[2], "Final Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
        header->ss[1], idate[0], month[idate[1] - 1], idate[2],
        idate[3], idate[4], idate[5] );
    for ( i = strlen( header->ttl[2] ); i < EPHCOM_MAXTTL; i++ )
        header->ttl[2][i] = ' ';
    header->ttl[0][EPHCOM_MAXTTL] = header->ttl[1][EPHCOM_MAXTTL] = header->ttl[2][EPHCOM_MAXTTL] = '\0';

//
//   ephcom_Header title lines.
//
//   Write the three title lines to the output file, padded with blanks,
//   84 characters long (84 is the first even multiple of 6 that is > 80,
//   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;
    }
//
//   Gregorian calendar.
//
    else // Explanatory Supplement to Astronomical Almanac, p. 604
    {
        L = j + 68569;
        N = ( 4 * L ) / 146097;
        L = L - ( 146097 * N + 3 ) / 4;
        I = ( 4000 * ( L + 1 ) ) / 1461001;
        L = L - ( 1461 * I ) / 4 + 31;
        J = ( 80 * L ) / 2447;
        D = L - ( 2447 * J ) / 80;
        L = J / 11;
        M = J + 2 - 12 * L;
        Y = 100 * ( N - 49 ) + I + L;
    }

    idate[0] = Y;
    idate[1] = M;
    idate[2] = D;
    idate[3] = ihour;
    idate[4] = imin;
    idate[5] = isec;
}

//! Convert calendar date and time to JD.  From pp. 604, 606 in the
//! Explanatory Supplement to the Astronomical Almanac.
//!
//! @param idate [IN ONLY] integer array of 6 values which contains
//! the integer year, month, day, hour, minute, and second.
//! @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.
//! @returns double-precision Julian Day number corresponding to idate.
//!
double
ephcom_cal2jd( int idate[6], int calendar_type )
{
    double tjd;
    int    jd;

//
//   Convert hours, minutes, seconds to fractional JD.
//
    tjd = ( idate[3] + ( idate[4] + idate[5] / 60.0 ) / 60.0 ) / 24.0 - 0.5;
//
//   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 &&
           ( idate[0] < 1582 ||                                // Before 1582
             ( idate[0] == 1582 &&
               ( idate[1] < 10 ||                              // Before October 1582
                 ( idate[1] == 10 && idate[2] < 15 ) ) ) ) ) ) // Before 15 October 1582
    {
        jd = 367 * idate[0] -
             ( 7 * ( idate[0] + 5001 + ( idate[1] - 9 ) / 7 ) ) / 4 +
             ( 275 * idate[1] ) / 9 +
             idate[2] + 1729777;
    }
//
//   Gregorian calendar.
//
    else // Explanatory Supplement to Astronomical Almanac, p. 604
    {
        jd = ( 1461 * ( idate[0] + 4800 + ( idate[1] - 14 ) / 12 ) ) / 4 +
             ( 367 * ( idate[1] - 2 - 12 * ( ( idate[1] - 14 ) / 12 ) ) ) / 12 -
             ( 3 * ( ( idate[0] + 4900 + ( idate[1] - 14 ) / 12 ) / 100 ) ) / 4 +
             idate[2] - 32075;
    }
//
//   Return value is whole JD number plus fractional JD number.
//
    tjd += (double) jd;

    return ( tjd );
}

//! If Julian date is close to an integer + 0.5, return that exact
//! "half" value.  Use of this routine makes the code more robust for
//! an early version of de422 (since corrected) which had numerical
//! noise in its half days.
//!
//! @param time [IN ONLY]Value of the Julian date which should be
//! exactly half integral.
//! @returns exactly half integral Julian date closest to time if
//! input time close to half-integral.  If input time is not close
//! to half-integral the routine returns the unmodified input time.
//!
static double
ephcom_exact_time( double time )
{
    double exact_half_time = ( time >= 0. ) ? (double) ( (int) time + 0.5 ) : (double) ( (int) time - 0.5 );
    if ( IF_SAME_DATE( time, exact_half_time ) )
        return exact_half_time;
    else
        return time;
}

//! Split time into an integer part and a positive remainder in the
//! semi-open range [0., 1.).
//!
//! @param time [IN ONLY]time value to be split.
//! @param integral_time [OUT ONLY]Pointer to a double value which
//! upon return will contain the integral part of the split time.
//! @returns remainder of the split time in the semi-open range [0., 1.).
//!

static double
ephcom_split( double time, double * integral_time )
{
    double retval = modf( time, integral_time );
    if ( retval < 0. )
    {
        retval         += 1.;
        *integral_time -= 1.;
    }
    return retval;
}