Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
//! @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.
//
Loading
Loading full blame…