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,