//! @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 and // for a given target and center, AFTER // calling ephcom_get_coords() (different // sequence than with JPL's FORTRAN PLEPH) // ephcom_get_coords() - calculate and // 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 #include #include #include #include // 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. //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //!
Index Identification
1 Mercury
2 Venus
3 Earth
4 Mars
5 Jupiter
6 Saturn
7 Uranus
8 Neptune
9 Pluto
10 Moon
11 Sun
12 Solar System Barycenter
13 Earth-Moon Barycenter
14 Nutation Angles
15 Libration Angles
16 Moon (Geocentric)
//! //! @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(): //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //!
Index Identification
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 Libration Angles
15 Moon (Geocentric)
//! //! @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: //! //! //! //! //! //! //! //! //! //! //! //! //! //! //! //!
Index Identification
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 Libration angles
//! //! @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; }