mirror of
https://github.com/cosinekitty/astronomy.git
synced 2025-12-30 11:09:33 -05:00
1039 lines
24 KiB
C
1039 lines
24 KiB
C
/*
|
|
Naval Observatory Vector Astrometry Software (NOVAS)
|
|
C Edition, Version 3.1
|
|
|
|
eph_manager.c: C version of JPL Ephemeris Manager for use with solsys1
|
|
|
|
U. S. Naval Observatory
|
|
Astronomical Applications Dept.
|
|
Washington, DC
|
|
http://www.usno.navy.mil/USNO/astronomical-applications
|
|
*/
|
|
|
|
#ifndef _EPHMAN_
|
|
#include "eph_manager.h"
|
|
#endif
|
|
|
|
/*
|
|
Define global variables
|
|
*/
|
|
|
|
short int KM;
|
|
|
|
/*
|
|
IPT and LPT defined as int to support 64 bit systems.
|
|
*/
|
|
|
|
int IPT[3][12], LPT[3];
|
|
|
|
long int NRL, NP, NV;
|
|
long int RECORD_LENGTH;
|
|
|
|
double SS[3], JPLAU, PC[18], VC[18], TWOT, EM_RATIO;
|
|
double *BUFFER;
|
|
|
|
FILE *EPHFILE = NULL;
|
|
|
|
/********ephem_open */
|
|
|
|
short int ephem_open (char *ephem_name,
|
|
|
|
double *jd_begin, double *jd_end,
|
|
short int *de_number)
|
|
/*
|
|
------------------------------------------------------------------------
|
|
|
|
PURPOSE:
|
|
This function opens a JPL planetary ephemeris file and
|
|
sets initial values. This function must be called
|
|
prior to calls to the other JPL ephemeris functions.
|
|
|
|
REFERENCES:
|
|
Standish, E.M. and Newhall, X X (1988). "The JPL Export
|
|
Planetary Ephemeris"; JPL document dated 17 June 1988.
|
|
|
|
INPUT
|
|
ARGUMENTS:
|
|
*ephem_name (char)
|
|
Name of the direct-access ephemeris file.
|
|
|
|
OUTPUT
|
|
ARGUMENTS:
|
|
*jd_begin (double)
|
|
Beginning Julian date of the ephemeris file.
|
|
*jd_end (double)
|
|
Ending Julian date of the ephemeris file.
|
|
*de_number (short int)
|
|
DE number of the ephemeris file opened.
|
|
|
|
RETURNED
|
|
VALUE:
|
|
(short int)
|
|
0 ...file exists and is opened correctly.
|
|
1 ...file does not exist/not found.
|
|
2-10...error reading from file header.
|
|
11 ...unable to set record length; ephemeris (DE number)
|
|
not in look-up table.
|
|
|
|
GLOBALS
|
|
USED:
|
|
SS eph_manager.h
|
|
JPLAU eph_manager.h
|
|
PC eph_manager.h
|
|
VC eph_manager.h
|
|
TWOT eph_manager.h
|
|
EM_RATIO eph_manager.h
|
|
BUFFER eph_manager.h
|
|
IPT eph_manager.h
|
|
LPT eph_manager.h
|
|
NRL eph_manager.h
|
|
KM eph_manager.h
|
|
NP eph_manager.h
|
|
NV eph_manager.h
|
|
RECORD_LENGTH eph_manager.h
|
|
EPHFILE eph_manager.h
|
|
|
|
FUNCTIONS
|
|
CALLED:
|
|
fclose stdio.h
|
|
free stdlib.h
|
|
fopen stdio.h
|
|
fread stdio.h
|
|
calloc stdlib.h
|
|
|
|
VER./DATE/
|
|
PROGRAMMER:
|
|
V1.0/06-90/JAB (USNO/NA)
|
|
V1.1/06-92/JAB (USNO/AA): Restructure and add initializations.
|
|
V1.2/07-98/WTH (USNO/AA): Modified to open files for different
|
|
ephemeris types. (200,403,404,405,406)
|
|
V1.3/11-07/WKP (USNO/AA): Updated prolog.
|
|
V1.4/09-10/WKP (USNO/AA): Changed ncon and denum variables and
|
|
sizeof ipt array to type 'int' for
|
|
64-bit system compatibility.
|
|
V1.5/09-10/WTH (USNO/AA): Added support for DE421, default case
|
|
for switch, close file on error.
|
|
V1.6/10-10/WKP (USNO/AA): Renamed function to lowercase to
|
|
comply with coding standards.
|
|
|
|
NOTES:
|
|
KM...flag defining physical units of the output states.
|
|
= 1, km and km/sec
|
|
= 0, AU and AU/day
|
|
Default value is 0 (KM determines time unit for nutations.
|
|
Angle unit is always radians.)
|
|
|
|
------------------------------------------------------------------------
|
|
*/
|
|
{
|
|
char ttl[252], cnam[2400];
|
|
|
|
short int i, j;
|
|
|
|
int ncon, denum;
|
|
|
|
ephem_close();
|
|
|
|
/*
|
|
Open file ephem_name.
|
|
*/
|
|
|
|
if ((EPHFILE = fopen (ephem_name, "rb")) == NULL)
|
|
{
|
|
return 1;
|
|
}
|
|
else
|
|
{
|
|
|
|
/*
|
|
File found. Set initializations and default values.
|
|
*/
|
|
|
|
KM = 0;
|
|
|
|
NRL = 0;
|
|
|
|
NP = 2;
|
|
NV = 3;
|
|
TWOT = 0.0;
|
|
|
|
for (i = 0; i < 18; i++)
|
|
{
|
|
PC[i] = 0.0;
|
|
VC[i] = 0.0;
|
|
}
|
|
|
|
PC[0] = 1.0;
|
|
VC[1] = 1.0;
|
|
|
|
/*
|
|
Read in values from the first record, aka the header.
|
|
*/
|
|
|
|
if (fread (ttl, sizeof ttl, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 2;
|
|
}
|
|
if (fread (cnam, sizeof cnam, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 3;
|
|
}
|
|
if (fread (SS, sizeof SS, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 4;
|
|
}
|
|
if (fread (&ncon, sizeof ncon, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 5;
|
|
}
|
|
if (fread (&JPLAU, sizeof JPLAU, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 6;
|
|
}
|
|
if (fread (&EM_RATIO, sizeof EM_RATIO, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 7;
|
|
}
|
|
for (i = 0; i < 12; i++)
|
|
for (j = 0; j < 3; j++)
|
|
if (fread (&IPT[j][i], sizeof(int), 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 8;
|
|
}
|
|
if (fread (&denum, sizeof denum, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 9;
|
|
}
|
|
if (fread (LPT, sizeof LPT, 1, EPHFILE) != 1)
|
|
{
|
|
ephem_close();
|
|
return 10;
|
|
}
|
|
|
|
/*
|
|
Set the value of the record length according to what JPL ephemeris is
|
|
being opened.
|
|
*/
|
|
|
|
switch (denum)
|
|
{
|
|
case 200:
|
|
RECORD_LENGTH = 6608;
|
|
break;
|
|
case 403: case 405:
|
|
case 421:
|
|
RECORD_LENGTH = 8144;
|
|
break;
|
|
case 404: case 406:
|
|
RECORD_LENGTH = 5824;
|
|
break;
|
|
|
|
/*
|
|
An unknown DE file was opened. Close the file and return an error
|
|
code.
|
|
*/
|
|
|
|
default:
|
|
*jd_begin = 0.0;
|
|
*jd_end = 0.0;
|
|
*de_number = 0;
|
|
ephem_close();
|
|
return 11;
|
|
}
|
|
|
|
BUFFER = (double *) calloc (RECORD_LENGTH / 8, sizeof(double));
|
|
|
|
*de_number = (short int) denum;
|
|
*jd_begin = SS[0];
|
|
*jd_end = SS[1];
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
/********ephem_close */
|
|
|
|
short int ephem_close (void)
|
|
/*
|
|
------------------------------------------------------------------------
|
|
|
|
PURPOSE:
|
|
This function closes a JPL planetary ephemeris file and frees the
|
|
memory.
|
|
|
|
REFERENCES:
|
|
Standish, E.M. and Newhall, X X (1988). "The JPL Export
|
|
Planetary Ephemeris"; JPL document dated 17 June 1988.
|
|
|
|
INPUT
|
|
ARGUMENTS:
|
|
None.
|
|
|
|
OUTPUT
|
|
ARGUMENTS:
|
|
None.
|
|
|
|
RETURNED
|
|
VALUE:
|
|
(short int)
|
|
0 ...file was already closed or closed correctly.
|
|
EOF...error closing file.
|
|
|
|
GLOBALS
|
|
USED:
|
|
BUFFER eph_manager.h
|
|
EPHFILE eph_manager.h
|
|
|
|
FUNCTIONS
|
|
CALLED:
|
|
fclose stdio.h
|
|
free stdlib.h
|
|
|
|
VER./DATE/
|
|
PROGRAMMER:
|
|
V1.0/11-07/WKP (USNO/AA)
|
|
V1.1/09-10/WKP (USNO/AA): Explicitly cast fclose return value to
|
|
type 'short int'.
|
|
V1.2/10-10/WKP (USNO/AA): Renamed function to lowercase to
|
|
comply with coding standards.
|
|
|
|
NOTES:
|
|
None.
|
|
|
|
------------------------------------------------------------------------
|
|
*/
|
|
{
|
|
short int error = 0;
|
|
|
|
if (EPHFILE)
|
|
{
|
|
error = (short int) fclose (EPHFILE);
|
|
EPHFILE = NULL; /* modified per errata in https://aa.usno.navy.mil/software/novas/novas_faq.php */
|
|
}
|
|
|
|
if (BUFFER)
|
|
{
|
|
free (BUFFER);
|
|
BUFFER = NULL;
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
/********planet_ephemeris */
|
|
|
|
short int planet_ephemeris (double tjd[2], short int target,
|
|
short int center,
|
|
|
|
double *position, double *velocity)
|
|
/*
|
|
------------------------------------------------------------------------
|
|
|
|
PURPOSE:
|
|
This function accesses the JPL planetary ephemeris to give the
|
|
position and velocity of the target object with respect to the
|
|
center object.
|
|
|
|
REFERENCES:
|
|
Standish, E.M. and Newhall, X X (1988). "The JPL Export
|
|
Planetary Ephemeris"; JPL document dated 17 June 1988.
|
|
|
|
INPUT
|
|
ARGUMENTS:
|
|
tjd[2] (double)
|
|
Two-element array containing the Julian date, which may be
|
|
split any way (although the first element is usually the
|
|
"integer" part, and the second element is the "fractional"
|
|
part). Julian date is in the TDB or "T_eph" time scale.
|
|
target (short int)
|
|
Number of 'target' point.
|
|
center (short int)
|
|
Number of 'center' (origin) point.
|
|
The numbering convention for 'target' and'center' is:
|
|
0 = Mercury 7 = Neptune
|
|
1 = Venus 8 = Pluto
|
|
2 = Earth 9 = Moon
|
|
3 = Mars 10 = Sun
|
|
4 = Jupiter 11 = Solar system bary.
|
|
5 = Saturn 12 = Earth-Moon bary.
|
|
6 = Uranus 13 = Nutations (long int. and obliq.)
|
|
(If nutations are desired, set 'target' = 13;
|
|
'center' will be ignored on that call.)
|
|
|
|
OUTPUT
|
|
ARGUMENTS:
|
|
*position (double)
|
|
Position vector array of target relative to center, measured
|
|
in AU.
|
|
*velocity (double)
|
|
Velocity vector array of target relative to center, measured
|
|
in AU/day.
|
|
|
|
RETURNED
|
|
VALUE:
|
|
(short int)
|
|
0 ...everything OK.
|
|
1,2...error returned from State.
|
|
|
|
GLOBALS
|
|
USED:
|
|
EM_RATIO eph_manager.h
|
|
|
|
FUNCTIONS
|
|
CALLED:
|
|
state eph_manager.h
|
|
|
|
VER./DATE/
|
|
PROGRAMMER:
|
|
V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
|
|
V1.1/07-93/WTH (USNO/AA): Update to C standards.
|
|
V2.0/07-98/WTH (USNO/AA): Modified for ease of use and linearity.
|
|
V3.0/11-06/JAB (USNO/AA): Allowed for use of input 'split' Julian
|
|
date for higher precision.
|
|
V3.1/11-07/WKP (USNO/AA): Updated prolog and error codes.
|
|
V3.1/12-07/WKP (USNO/AA): Removed unreferenced variables.
|
|
V3.2/10-10/WKP (USNO/AA): Renamed function to lowercase to
|
|
comply with coding standards.
|
|
|
|
NOTES:
|
|
None.
|
|
|
|
------------------------------------------------------------------------
|
|
*/
|
|
{
|
|
short int i, error = 0, earth = 2, moon = 9;
|
|
short int do_earth = 0, do_moon = 0;
|
|
|
|
double jed[2];
|
|
double pos_moon[3] = {0.0,0.0,0.0}, vel_moon[3] = {0.0,0.0,0.0},
|
|
pos_earth[3] = {0.0,0.0,0.0}, vel_earth[3] = {0.0,0.0,0.0};
|
|
double target_pos[3] = {0.0,0.0,0.0}, target_vel[3] = {0.0,0.0,0.0},
|
|
center_pos[3] = {0.0,0.0,0.0}, center_vel[3] = {0.0,0.0,0.0};
|
|
|
|
/*
|
|
Initialize 'jed' for 'state' and set up component count.
|
|
*/
|
|
|
|
jed[0] = tjd[0];
|
|
jed[1] = tjd[1];
|
|
|
|
/*
|
|
Check for target point = center point.
|
|
*/
|
|
|
|
if (target == center)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
position[i] = 0.0;
|
|
velocity[i] = 0.0;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
/*
|
|
Check for instances of target or center being Earth or Moon,
|
|
and for target or center being the Earth-Moon barycenter.
|
|
*/
|
|
|
|
if ((target == earth) || (center == earth))
|
|
do_moon = 1;
|
|
if ((target == moon) || (center == moon))
|
|
do_earth = 1;
|
|
if ((target == 12) || (center == 12))
|
|
do_earth = 1;
|
|
|
|
if (do_earth)
|
|
{
|
|
error = state (jed,2, pos_earth,vel_earth);
|
|
if (error)
|
|
return error;
|
|
}
|
|
|
|
if (do_moon)
|
|
{
|
|
error = state (jed,9, pos_moon,vel_moon);
|
|
if (error)
|
|
return error;
|
|
}
|
|
|
|
/*
|
|
Make call to State for target object.
|
|
*/
|
|
|
|
if (target == 11)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
target_pos[i] = 0.0;
|
|
target_vel[i] = 0.0;
|
|
}
|
|
}
|
|
else if (target == 12)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
target_pos[i] = pos_earth[i];
|
|
target_vel[i] = vel_earth[i];
|
|
}
|
|
}
|
|
else
|
|
error = state (jed,target, target_pos,target_vel);
|
|
|
|
if (error)
|
|
return error;
|
|
|
|
/*
|
|
Make call to State for center object.
|
|
*/
|
|
|
|
/*
|
|
If the requested center is the Solar System barycenter,
|
|
then don't bother with a second call to State.
|
|
*/
|
|
|
|
if (center == 11)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
center_pos[i] = 0.0;
|
|
center_vel[i] = 0.0;
|
|
}
|
|
}
|
|
|
|
/*
|
|
Center is Earth-Moon barycenter, which was already computed above.
|
|
*/
|
|
|
|
else if (center == 12)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
center_pos[i] = pos_earth[i];
|
|
center_vel[i] = vel_earth[i];
|
|
}
|
|
}
|
|
else
|
|
error = state (jed,center, center_pos,center_vel);
|
|
|
|
if (error)
|
|
return error;
|
|
|
|
/*
|
|
Check for cases of Earth as target and Moon as center or vice versa.
|
|
*/
|
|
|
|
if ((target == earth) && (center == moon))
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
position[i] = -center_pos[i];
|
|
velocity[i] = -center_vel[i];
|
|
}
|
|
return 0;
|
|
}
|
|
else if ((target == moon) && (center == earth))
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
position[i] = target_pos[i];
|
|
velocity[i] = target_vel[i];
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
/*
|
|
Check for Earth as target, or as center.
|
|
*/
|
|
|
|
else if (target == earth)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
target_pos[i] = target_pos[i] - (pos_moon[i] /
|
|
(1.0 + EM_RATIO));
|
|
target_vel[i] = target_vel[i] - (vel_moon[i] /
|
|
(1.0 + EM_RATIO));
|
|
}
|
|
}
|
|
else if (center == earth)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
center_pos[i] = center_pos[i] - (pos_moon[i] /
|
|
(1.0 + EM_RATIO));
|
|
center_vel[i] = center_vel[i] - (vel_moon[i] /
|
|
(1.0 + EM_RATIO));
|
|
}
|
|
}
|
|
|
|
/*
|
|
Check for Moon as target, or as center.
|
|
*/
|
|
|
|
else if (target == moon)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
target_pos[i] = (pos_earth[i] - (target_pos[i] /
|
|
(1.0 + EM_RATIO))) + target_pos[i];
|
|
target_vel[i] = (vel_earth[i] - (target_vel[i] /
|
|
(1.0 + EM_RATIO))) + target_vel[i];
|
|
}
|
|
}
|
|
else if (center == moon)
|
|
{
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
center_pos[i] = (pos_earth[i] - (center_pos[i] /
|
|
(1.0 + EM_RATIO))) + center_pos[i];
|
|
center_vel[i] = (vel_earth[i] - (center_vel[i] /
|
|
(1.0 + EM_RATIO))) + center_vel[i];
|
|
}
|
|
}
|
|
|
|
/*
|
|
Compute position and velocity vectors.
|
|
*/
|
|
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
position[i] = target_pos[i] - center_pos[i];
|
|
velocity[i] = target_vel[i] - center_vel[i];
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
/********state */
|
|
|
|
short int state (double *jed, short int target,
|
|
|
|
double *target_pos, double *target_vel)
|
|
/*
|
|
------------------------------------------------------------------------
|
|
|
|
PURPOSE:
|
|
This function reads and interpolates the JPL planetary
|
|
ephemeris file.
|
|
|
|
REFERENCES:
|
|
Standish, E.M. and Newhall, X X (1988). "The JPL Export
|
|
Planetary Ephemeris"; JPL document dated 17 June 1988.
|
|
|
|
INPUT
|
|
ARGUMENTS:
|
|
*jed (double)
|
|
2-element Julian date (TDB) at which interpolation is wanted.
|
|
Any combination of jed[0]+jed[1] which falls within the time
|
|
span on the file is a permissible epoch. See Note 1 below.
|
|
target (short int)
|
|
The requested body to get data for from the ephemeris file.
|
|
The designation of the astronomical bodies is:
|
|
0 = Mercury 6 = Uranus
|
|
1 = Venus 7 = Neptune
|
|
2 = Earth-Moon barycenter 8 = Pluto
|
|
3 = Mars 9 = geocentric Moon
|
|
4 = Jupiter 10 = Sun
|
|
5 = Saturn
|
|
|
|
OUTPUT
|
|
ARGUMENTS:
|
|
*target_pos (double)
|
|
The barycentric position vector array of the requested object,
|
|
in AU.
|
|
(If target object is the Moon, then the vector is geocentric.)
|
|
*target_vel (double)
|
|
The barycentric velocity vector array of the requested object,
|
|
in AU/Day.
|
|
|
|
Both vectors are referenced to the Earth mean equator and
|
|
equinox of epoch.
|
|
|
|
RETURNED
|
|
VALUE:
|
|
(short int)
|
|
0...everything OK.
|
|
1...error reading ephemeris file.
|
|
2...epoch out of range.
|
|
3...target out of range.
|
|
|
|
GLOBALS
|
|
USED:
|
|
KM eph_manager.h
|
|
EPHFILE eph_manager.h
|
|
IPT eph_manager.h
|
|
BUFFER eph_manager.h
|
|
NRL eph_manager.h
|
|
RECORD_LENGTH eph_manager.h
|
|
SS eph_manager.h
|
|
JPLAU eph_manager.h
|
|
|
|
FUNCTIONS
|
|
CALLED:
|
|
split eph_manager.h
|
|
fseek stdio.h
|
|
fread stdio.h
|
|
interpolate eph_manager.h
|
|
ephem_close eph_manager.h
|
|
|
|
VER./DATE/
|
|
PROGRAMMER:
|
|
V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
|
|
V1.1/07-93/WTH (USNO/AA): Update to C standards.
|
|
V2.0/07-98/WTH (USNO/AA): Modify to make position and velocity
|
|
two distinct vector arrays. Routine set
|
|
to compute one state per call.
|
|
V2.1/11-07/WKP (USNO/AA): Updated prolog.
|
|
V2.2/10-10/WKP (USNO/AA): Renamed function to lowercase to
|
|
comply with coding standards.
|
|
|
|
NOTES:
|
|
1. For ease in programming, the user may put the entire epoch in
|
|
jed[0] and set jed[1] = 0. For maximum interpolation accuracy,
|
|
set jed[0] = the most recent midnight at or before
|
|
interpolation epoch, and set jed[1] = fractional part of a day
|
|
elapsed between jed[0] and epoch. As an alternative, it may
|
|
prove convenient to set jed[0] = some fixed epoch, such as
|
|
start of the integration and jed[1] = elapsed interval between
|
|
then and epoch.
|
|
|
|
------------------------------------------------------------------------
|
|
*/
|
|
{
|
|
short int i;
|
|
|
|
long int nr, rec;
|
|
|
|
double t[2], aufac = 1.0, jd[4], s;
|
|
|
|
/*
|
|
Check for invalid targets.
|
|
*/
|
|
if (target < 0 || target > 10)
|
|
return 3;
|
|
|
|
/*
|
|
Set units based on value of the 'KM' flag.
|
|
*/
|
|
|
|
if (KM)
|
|
t[1] = SS[2] * 86400.0;
|
|
else
|
|
{
|
|
t[1] = SS[2];
|
|
aufac = 1.0 / JPLAU;
|
|
}
|
|
|
|
/*
|
|
Check epoch.
|
|
*/
|
|
|
|
s = jed[0] - 0.5;
|
|
split (s, &jd[0]);
|
|
split (jed[1], &jd[2]);
|
|
jd[0] += jd[2] + 0.5;
|
|
jd[1] += jd[3];
|
|
split (jd[1], &jd[2]);
|
|
jd[0] += jd[2];
|
|
|
|
/*
|
|
Return error code if date is out of range.
|
|
*/
|
|
|
|
if ((jd[0] < SS[0]) || ((jd[0] + jd[3]) > SS[1]))
|
|
return 2;
|
|
|
|
/*
|
|
Calculate record number and relative time interval.
|
|
*/
|
|
|
|
nr = (long int) ((jd[0] - SS[0]) / SS[2]) + 3;
|
|
if (jd[0] == SS[1])
|
|
nr -= 2;
|
|
t[0] = ((jd[0] - ((double) (nr-3) * SS[2] + SS[0])) + jd[3]) / SS[2];
|
|
|
|
/*
|
|
Read correct record if it is not already in memory.
|
|
*/
|
|
|
|
if (nr != NRL)
|
|
{
|
|
NRL = nr;
|
|
rec = (nr - 1) * RECORD_LENGTH;
|
|
fseek (EPHFILE, rec, SEEK_SET);
|
|
if (!fread (BUFFER, RECORD_LENGTH, 1, EPHFILE))
|
|
{
|
|
ephem_close ();
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
/*
|
|
Check and interpolate for requested body.
|
|
*/
|
|
|
|
interpolate (&BUFFER[IPT[0][target]-1],t,IPT[1][target],
|
|
IPT[2][target], target_pos,target_vel);
|
|
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
target_pos[i] *= aufac;
|
|
target_vel[i] *= aufac;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
/********interpolate */
|
|
|
|
void interpolate (double *buf, double *t, long int ncf, long int na,
|
|
|
|
double *position, double *velocity)
|
|
/*
|
|
------------------------------------------------------------------------
|
|
|
|
PURPOSE:
|
|
This function differentiates and interpolates a set of
|
|
Chebyshev coefficients to give position and velocity.
|
|
|
|
REFERENCES:
|
|
Standish, E.M. and Newhall, X X (1988). "The JPL Export
|
|
Planetary Ephemeris"; JPL document dated 17 June 1988.
|
|
|
|
INPUT
|
|
ARGUMENTS:
|
|
*buf (double)
|
|
Array of Chebyshev coefficients of position.
|
|
*t (double)
|
|
t[0] is fractional time interval covered by coefficients at
|
|
which interpolation is desired (0 <= t[0] <= 1).
|
|
t[1] is length of whole interval in input time units.
|
|
ncf (long int)
|
|
Number of coefficients per component.
|
|
na (long int)
|
|
Number of sets of coefficients in full array
|
|
(i.e., number of sub-intervals in full interval).
|
|
|
|
OUTPUT
|
|
ARGUMENTS:
|
|
*position (double)
|
|
Position array of requested object.
|
|
*velocity (double)
|
|
Velocity array of requested object.
|
|
|
|
RETURNED
|
|
VALUE:
|
|
None.
|
|
|
|
GLOBALS
|
|
USED:
|
|
NP eph_manager.h
|
|
NV eph_manager.h
|
|
PC eph_manager.h
|
|
VC eph_manager.h
|
|
TWOT eph_manager.h
|
|
|
|
FUNCTIONS
|
|
CALLED:
|
|
fmod math.h
|
|
|
|
VER./DATE/
|
|
PROGRAMMER:
|
|
V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
|
|
V1.1/07-93/WTH (USNO/AA): Update to C standards.
|
|
V1.2/07-98/WTH (USNO/AA): Modify to make position and velocity
|
|
two distinct vector arrays.
|
|
V1.3/11-07/WKP (USNO/AA): Updated prolog.
|
|
V1.4/12-07/WKP (USNO/AA): Changed ncf and na arguments from short
|
|
int to long int.
|
|
V1.5/10-10/WKP (USNO/AA): Renamed function to lowercase to
|
|
comply with coding standards.
|
|
|
|
NOTES:
|
|
None.
|
|
|
|
------------------------------------------------------------------------
|
|
*/
|
|
{
|
|
long int i, j, k, l;
|
|
|
|
double dna, dt1, temp, tc, vfac;
|
|
|
|
/*
|
|
Get correct sub-interval number for this set of coefficients and
|
|
then get normalized Chebyshev time within that subinterval.
|
|
*/
|
|
|
|
dna = (double) na;
|
|
dt1 = (double) ((long int) t[0]);
|
|
temp = dna * t[0];
|
|
l = (long int) (temp - dt1);
|
|
|
|
/*
|
|
'tc' is the normalized Chebyshev time (-1 <= tc <= 1).
|
|
*/
|
|
|
|
tc = 2.0 * (fmod (temp, 1.0) + dt1) - 1.0;
|
|
|
|
/*
|
|
Check to see whether Chebyshev time has changed, and compute new
|
|
polynomial values if it has. (The element PC[1] is the value of
|
|
t1[tc] and hence contains the value of 'tc' on the previous call.)
|
|
*/
|
|
|
|
if (tc != PC[1])
|
|
{
|
|
NP = 2;
|
|
NV = 3;
|
|
PC[1] = tc;
|
|
TWOT = tc + tc;
|
|
}
|
|
|
|
/*
|
|
Be sure that at least 'ncf' polynomials have been evaluated and
|
|
are stored in the array 'PC'.
|
|
*/
|
|
|
|
if (NP < ncf)
|
|
{
|
|
for (i = NP; i < ncf; i++)
|
|
PC[i] = TWOT * PC[i-1] - PC[i-2];
|
|
NP = ncf;
|
|
}
|
|
|
|
/*
|
|
Interpolate to get position for each component.
|
|
*/
|
|
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
position[i] = 0.0;
|
|
for (j = ncf-1; j >= 0; j--)
|
|
{
|
|
k = j + (i * ncf) + (l * (3 * ncf));
|
|
position[i] += PC[j] * buf[k];
|
|
}
|
|
}
|
|
|
|
/*
|
|
If velocity interpolation is desired, be sure enough derivative
|
|
polynomials have been generated and stored.
|
|
*/
|
|
|
|
vfac = (2.0 * dna) / t[1];
|
|
VC[2] = 2.0 * TWOT;
|
|
if (NV < ncf)
|
|
{
|
|
for (i = NV; i < ncf; i++)
|
|
VC[i] = TWOT * VC[i-1] + PC[i-1] + PC[i-1] - VC[i-2];
|
|
NV = ncf;
|
|
}
|
|
|
|
/*
|
|
Interpolate to get velocity for each component.
|
|
*/
|
|
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
velocity[i] = 0.0;
|
|
for (j = ncf-1; j > 0; j--)
|
|
{
|
|
k = j + (i * ncf) + (l * (3 * ncf));
|
|
velocity[i] += VC[j] * buf[k];
|
|
}
|
|
velocity[i] *= vfac;
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
/********split */
|
|
|
|
void split (double tt,
|
|
|
|
double *fr)
|
|
/*
|
|
------------------------------------------------------------------------
|
|
|
|
PURPOSE:
|
|
This function breaks up a double number into a double integer
|
|
part and a fractional part.
|
|
|
|
REFERENCES:
|
|
Standish, E.M. and Newhall, X X (1988). "The JPL Export
|
|
Planetary Ephemeris"; JPL document dated 17 June 1988.
|
|
|
|
INPUT
|
|
ARGUMENTS:
|
|
tt (double)
|
|
Input number.
|
|
|
|
OUTPUT
|
|
ARGUMENTS:
|
|
*fr (double)
|
|
2-element output array;
|
|
fr[0] contains integer part,
|
|
fr[1] contains fractional part.
|
|
For negative input numbers,
|
|
fr[0] contains the next more negative integer;
|
|
fr[1] contains a positive fraction.
|
|
|
|
RETURNED
|
|
VALUE:
|
|
None.
|
|
|
|
GLOBALS
|
|
USED:
|
|
None.
|
|
|
|
FUNCTIONS
|
|
CALLED:
|
|
None.
|
|
|
|
VER./DATE/
|
|
PROGRAMMER:
|
|
V1.0/06-90/JAB (USNO/NA): CA coding standards
|
|
V1.1/03-93/WTH (USNO/AA): Convert to C.
|
|
V1.2/07-93/WTH (USNO/AA): Update to C standards.
|
|
V1.3/10-10/WKP (USNO/AA): Renamed function to lowercase to
|
|
comply with coding standards.
|
|
|
|
NOTES:
|
|
None.
|
|
|
|
------------------------------------------------------------------------
|
|
*/
|
|
{
|
|
|
|
/*
|
|
Get integer and fractional parts.
|
|
*/
|
|
|
|
fr[0] = (double)((long int) tt);
|
|
fr[1] = tt - fr[0];
|
|
|
|
/*
|
|
Make adjustments for negative input number.
|
|
*/
|
|
|
|
if ((tt >= 0.0) || (fr[1] == 0.0))
|
|
return;
|
|
else
|
|
{
|
|
fr[0] = fr[0] - 1.0;
|
|
fr[1] = fr[1] + 1.0;
|
|
}
|
|
|
|
return;
|
|
}
|