ephcom.h

//! @file
//! Header information 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.

#ifndef __EPHCOM_H__
#define __EPHCOM_H__

// Set up for dealing with function visibility issues.
#include "ephcomdll.h"

#define EPHCOM_VERSION       "1.0"

#define EPHCOM_MAXLINE       128 // Maximum # characters to allow in input line
#define EPHCOM_MAXTTL        84  // Maximum # characters of information in ttl[0-3]
#define EPHCOM_MINJD         -999999999.5
#define EPHCOM_MAXJD         999999999.5

//   Objects for pleph() ntarget and ncenter parameters.
#define EPHCOM_MERCURY       1
#define EPHCOM_VENUS         2
#define EPHCOM_EARTH         3
#define EPHCOM_MARS          4
#define EPHCOM_JUPITER       5
#define EPHCOM_SATURN        6
#define EPHCOM_URANUS        7
#define EPHCOM_NEPTUNE       8
#define EPHCOM_PLUTO         9
#define EPHCOM_MOON          10    // Moon, relative to Solar System center
#define EPHCOM_SUN           11
#define EPHCOM_SSBARY        12    // Solar System Barycenter
#define EPHCOM_EMBARY        13    // Earth-Moon Barycenter
#define EPHCOM_NUTATION      14
#define EPHCOM_LIBRATION     15
#define EPHCOM_GEOMOON       16    // Original Lunar Ephemeris coordinates
#define EPHCOM_OBSERVER      17    // User-defined observer position
#define EPHCOM_NUMOBJECTS    17    // Allocate memory for 17 solar sys objs

//! This structure holds all the information contained in a JPLEPH header.
//! When an ASCII or binary header is read, this structure is populated.
//! Fill out this structure before writing an ASCII or binary header, and
//! before performing any interpolations.
//!
#ifdef __cplusplus
extern "C"{
    #endif
typedef struct
{
    //! block size, in first line of ASCII header.
    int    ksize;
    //! Number of Chebyshev coefficients in data blocks.
    int    ncoeff;
    //! Hold up to 14*6=EPHCOM_MAXTTL characters + "\n" + null.
    char   ttl[3][EPHCOM_MAXTTL + 2];
    //! Number of assigned in cnam.
    int    ncon;
    //! Hold up to 400 6-character names ending with null.
    char   cnam[400][7];
    //! Number of values for cval, to compare with ncon.
    int    nval;
    //! Constant values, corresponding to cnam names.
    double cval[400];
    //! Astronomical unit in km.
    double au;
    //! Earth-Moon mass ratio.
    double emrat;
    //! Speed of light, km/sec.
    double clight;
    //! Ephemeris number.
    int    numde;
    //! Lunar ephemeris number (can be same # as numde).
    int    numle;
    //! Start Julian day number, stop Julian day number, and step size (in days).
    double ss[3];
    //! Index pointers into Chebyshev coefficients.
    int    ipt[12][3];
    //! Libration pointer in a block.
    int    lpt[3];
    //! Maximum Chebyshev coefficients for a body.
    int    maxcheby;
} ephcom_Header;

//!   This struct holds all interpolated positions of planets, Sun, and Moon
//!   at a given time.  All of the information available from interpolation
//!   at a given point in time is computed at once and preserved here.
//!
//!   To populate this structure, have the ephemeris file open, and set:
//!
//!   km, seconds, bary, et2[0], and et2[1] in the ephcom_Coords struct whose
//!   pointer is an argument to ephcom_get_coords().
//!
//!   Then call ephcom_get_coords() to get all coordinates, then call
//!   ephcom_pleph() for each desired (ntarget,ncenter) combination.
//!   See testeph.c for an example.  Note that unlike JPL's PLEPH()
//!   subroutine, you cannot call ephcom_pleph() without first initializing
//!   the pv[] array in ephcom_get_coords().
//!
//!   The first index of pv is
//!   interpreted as follows (note the offset by one compared to the
//!   interpretation of the index arguments to 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>
//!   <tr> <td>16</td> <td>Possible user-defined coordinates</td> </tr>
//! </table>
//!
//!   There are some extra entries at the end, compared to JPL's PLEPH():
//!
//!      pv[15][] - preserves the unmodified lunar ephemeris coordinates.
//!                 These coordinates are never offset by center, so that
//!                 their full precision is maintained.
//!      pv[16][] - start of possible user-defined coordinates.
//!                 These coordinates are offset by center.
//!
typedef struct
{
    //! 1 = positions in km; 0 = positions in AU.
    int    km;
    //! 1 = timescale is seconds; 0 = timescale is days.
    int    seconds;
    //! 1 = Barycentric coordinates; 0 = adjust for center.
    int    bary;
    //! object to use as center (instead of SSBARY).
    int    center;
    //! 0 = raw JPL ephemeris coordinates. Note: none of the supplied routines read or modify the coordtype value.
    int    coordtype;
    //! Julian Day of interpolation. For best precision, et2[0] should be an exact integral or exact half-integral number of days while et2[1] should be a correction to et2[0] between 0. and 1.
    double et2[2];
    //! x, y, z Position & Velocity.  The first index of pv is interpreted as in the above table.
    double pv[EPHCOM_NUMOBJECTS][6];
} ephcom_Coords;

// Useful macro for Julian date comparison.
// JPL ephemeris Julian dates (so far) range from 6.e5 (3000 BC) to 2.8e6 (3000 AD).  So
// double precision (64-bit floating point) should be able to represent a date to within
// ~ 1.e-15*2.8e6 = 2.8e-9 days (or 0.24 ms).  Make criterion ~100 times that value since
// all JPL ephemerides up to now have dates which are offset by 0.5 from an integer so
// that a criterion of 0.25 days would probably work.

#define JULIAN_DATE_CRITERION    ( 3.e-7 )

// Need this logic to compare dates rather than exact equalities because
// (for DE422 at least) rounding errors have crept into the dates.
#define IF_SAME_DATE( a, b )    ( fabs( a - b ) <= JULIAN_DATE_CRITERION )

// Here is where the public API for libephcom is described.

// Read a JPL Ephemeris ASCII header from the FILE pointed to by the
// infp argument and store and return values in the ephcom_Header
// struct pointed to by the header argument.  Write any errors to
// stderr.
EPHCOMDLLIMPEXP void
ephcom_readascii_header( FILE * infp, ephcom_Header *header );

//
//   Read a block of data coefficients from a JPL ASCII ephemeris file.
//   Returns number of coefficients read, 0 at EOF.
//
EPHCOMDLLIMPEXP int
ephcom_readascii_block(
    FILE * infp,
    ephcom_Header *header,
    double *datablock );

//
//   Read a JPL Ephemeris header in binary format.  Store values in
//   an ephcom_Header struct.
//
EPHCOMDLLIMPEXP void
ephcom_readbinary_header( FILE * infp, ephcom_Header *header );

//
//   Read a JPL Ephemeris data block in binary format.
//
//   This is the only routine in this library that accesses a file
//   as a direct access file, with a specified block number.  The
//   block number ranges from 0 on up (starting at first data block,
//   after the 2 header blocks).  Returns the number of coefficients
//   read, or 0 at EOF.
//
EPHCOMDLLIMPEXP int
ephcom_readbinary_block(
    FILE *infp,            // File pointer for direct access file
    ephcom_Header *header, // header struct, already filled in
    int blocknum,          // Data block number, starting with 0
    double *datablock      // returned coefficient data block
    );

//
//   Write header information in ASCII format.
//
EPHCOMDLLIMPEXP void
ephcom_writeascii_header( FILE * outfp, ephcom_Header *header );

//
//   Write coefficient block information in ASCII format.
//
EPHCOMDLLIMPEXP void
ephcom_writeascii_block(
    FILE * outfp,
    ephcom_Header *header,
    int blocknum,
    double *datablock );

//
//   Write a JPL Ephemeris header in binary format.
//
EPHCOMDLLIMPEXP void
ephcom_writebinary_header( FILE *outfp, ephcom_Header *header );

//
//   Write a block of data coefficients in JPL binary file format.
//
EPHCOMDLLIMPEXP void
ephcom_writebinary_block(
    FILE * outfp,
    ephcom_Header *header,
    int blocknum,
    double *datablock );

//
//   ephcom_parse_block() - Parse a binary block of data.  Warning: verbose!
//                          Writes parsed output to file pointer outfp.
//
EPHCOMDLLIMPEXP void
ephcom_parse_block(
    FILE * outfp,
    ephcom_Header *header,
    double *datablock );

//
//   Read in a double precision value from the given file with bytes swapped
//   if necessary to match host order (Little- or DEC- Endian).  On Intel 80x86
//   the bytes will get swapped, on Motorola or SPARC they won't.
//
EPHCOMDLLIMPEXP double
ephcom_indouble( FILE *infp );

//
//   Read in an integer (4--byte) value to the given file with bytes swapped
//   if necessary to match host order (Little- or DEC- Endian).  On Intel 80x86
//   the bytes will get swapped, on Motorola or SPARC they won't.
//
EPHCOMDLLIMPEXP int
ephcom_inint( FILE *infp );

//
//   Planetary Ephemeris.  Takes coordinates already calculated in
//   coords structure and converts to vectors and vector dot in testr[].
//   Bodies start at 1 for Mercury, to match the JPL PLEPH() numbering.
//   Values for ntarg and ncntr correspond to locations ntarg-1 and
//   ncntr-1 in coords->pv[].
//
EPHCOMDLLIMPEXP void
ephcom_pleph( ephcom_Coords *coords, int ntarg, int ncntr, double *r );

//
//   ephcom_get_coords() - Interpolate positions and velocities at given time.
//
EPHCOMDLLIMPEXP int
ephcom_get_coords( FILE * infp,
                   ephcom_Header *header,
                   ephcom_Coords *coords,
                   double *datablock );

//
//   ephcom_cal2jd() - convert calendar date and time to JD.
//
//      idate: integer year, month, day, hour, minute, second
//      calendar_type: -1=Julian; 0=Automatic; 1=Gregorian
//      return value: double precision Julian Day of idate[]
//
//   From pp. 604, 606 in the Explanatory Supplement to the Astronomical Almanac.
//
EPHCOMDLLIMPEXP double
ephcom_cal2jd( int idate[6], int calendar_type );
#ifdef __cplusplus
}
#endif
#endif  // __EPHCOM_H__