mirror
/
XEphem
mirror of https://github.com/XEphem/XEphem.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
XEphem/libastro/vsop87.c

208 lines
6.0 KiB
C
Raw Permalink Blame History

/* VSOP87 planetary theory
*
* currently uses version VSOP87D:
* heliocentric spherical, mean ecliptic of date.
*
* calculation of rates (daily changes) is optional;
* see header file for the necessary #define's
*
* rough orientation on calculation time, miliseconds
* on an HP 715/75, all planets Mercury to Neptune, prec=0.0:
*
* terms with rates without rates
* 3598 11 7.1
* 31577 51 44
*
* with secular terms for JD 2232395.0 19/12/1399 0h TDB:
*
* FULL PRECISION code (31577 terms), milliseconds
* prec terms rates no rates
* 1e-8 15086 62 36
* 1e-7 10105 44 25
* 1e-6 3725 20 13
* 1e-5 1324 11 7.8
* 1e-4 443 7.0 6.0
* 1e-3 139 6.0 5.0
*
* REDUCED PRECISION code (3598 terms), milliseconds
* prec terms rates no rates
* 1e-7 2463 9.9 5.5
* 1e-6 1939 8.0 4.5
* 1e-5 1131 4.9 2.9
* 1e-4 443 2.2 1.5
* 1e-3 139 1.0 0.9
*/
#include <math.h>
#include "astro.h"
#include "vsop87.h"
#define VSOP_A1000 365250.0 /* days per millenium */
#define VSOP_MAXALPHA 5 /* max degree of time */
/******************************************************************
* adapted from BdL FORTRAN Code; stern
*
* Reference : Bureau des Longitudes - PBGF9502
*
* Object : calculate a VSOP87 position for a given time.
*
* Input :
*
* mj modified julian date, counted from J1900.0
* time scale : dynamical time TDB.
*
* obj object number as in astro.h, NB: not for pluto
*
* prec relative precision
*
* if prec is equal to 0 then the precision is the precision
* p0 of the complete solution VSOP87.
* Mercury p0 = 0.6 10**-8
* Venus p0 = 2.5 10**-8
* Earth p0 = 2.5 10**-8
* Mars p0 = 10.0 10**-8
* Jupiter p0 = 35.0 10**-8
* Saturn p0 = 70.0 10**-8
* Uranus p0 = 8.0 10**-8
* Neptune p0 = 42.0 10**-8
*
* if prec is not equal to 0, let us say in between p0 and
* 10**-3, the precision is :
* for the positions :
* - prec*a0 au for the distances.
* - prec rad for the other variables.
* for the velocities :
* - prec*a0 au/day for the distances.
* - prec rad/day for the other variables.
* a0 is the semi-major axis of the body.
*
* Output :
*
* ret[6] array of the results (double).
*
* for spherical coordinates :
* 1: longitude (rd)
* 2: latitude (rd)
* 3: radius (au)
* #if VSOP_GETRATE:
* 4: longitude velocity (rad/day)
* 5: latitude velocity (rad/day)
* 6: radius velocity (au/day)
*
* return: error index (int)
* 0: no error.
* 2: object out of range [MERCURY .. NEPTUNE, SUN]
* 3: precision out of range [0.0 .. 1e-3]
******************************************************************/
int
vsop87 (double mj, int obj, double prec, double *ret)
{
static double (*vx_map[])[3] = { /* data tables */
vx_mercury, vx_venus, vx_mars, vx_jupiter,
vx_saturn, vx_uranus, vx_neptune, 0, vx_earth,
};
static int (*vn_map[])[3] = { /* indexes */
vn_mercury, vn_venus, vn_mars, vn_jupiter,
vn_saturn, vn_uranus, vn_neptune, 0, vn_earth,
};
static double a0[] = { /* semimajor axes; for precision ctrl only */
0.39, 0.72, 1.5, 5.2, 9.6, 19.2, 30.1, 39.5, 1.0,
};
double (*vx_obj)[3] = vx_map[obj]; /* VSOP87 data and indexes */
int (*vn_obj)[3] = vn_map[obj];
double t[VSOP_MAXALPHA+1]; /* powers of time */
double t_abs[VSOP_MAXALPHA+1]; /* powers of abs(time) */
double q; /* aux for precision control */
int i, cooidx, alpha; /* misc indexes */
if (obj == PLUTO || obj > SUN)
return (2);
if (prec < 0.0 || prec > 1e-3)
return(3);
/* zero result array */
for (i = 0; i < 6; ++i) ret[i] = 0.0;
/* time and its powers */
t[0] = 1.0;
t[1] = (mj - J2000)/VSOP_A1000;
for (i = 2; i <= VSOP_MAXALPHA; ++i) t[i] = t[i-1] * t[1];
t_abs[0] = 1.0;
for (i = 1; i <= VSOP_MAXALPHA; ++i) t_abs[i] = fabs(t[i]);
/* precision control */
q = -log10(prec + 1e-35) - 2; /* decades below 1e-2 */
q = VSOP_ASCALE * prec / 10.0 / q; /* reduce threshold progressively
* for higher precision */
/* do the term summation; first the spatial dimensions */
for (cooidx = 0; cooidx < 3; ++cooidx) {
/* then the powers of time */
for (alpha = 0; vn_obj[alpha+1][cooidx] ; ++alpha) {
double p, term, termdot;
/* precision threshold */
p= alpha ? q/(t_abs[alpha] + alpha*t_abs[alpha-1]*1e-4 + 1e-35) : q;
#if VSOP_SPHERICAL
if (cooidx == 2) /* scale by semimajor axis for radius */
#endif
p *= a0[obj];
term = termdot = 0.0;
for (i = vn_obj[alpha][cooidx]; i < vn_obj[alpha+1][cooidx]; ++i) {
double a, b, c, arg;
a = vx_obj[i][0];
if (a < p) continue; /* ignore small terms */
b = vx_obj[i][1];
c = vx_obj[i][2];
arg = b + c * t[1];
term += a * cos(arg);
#if VSOP_GETRATE
termdot += -c * a * sin(arg);
#endif
}
ret[cooidx] += t[alpha] * term;
#if VSOP_GETRATE
ret[cooidx + 3] += t[alpha] * termdot +
((alpha > 0) ? alpha * t[alpha - 1] * term : 0.0);
#endif
} /* alpha */
} /* cooidx */
for (i = 0; i < 6; ++i) ret[i] /= VSOP_ASCALE;
#if VSOP_SPHERICAL
/* reduce longitude to 0..2pi */
ret[0] -= floor(ret[0]/(2.*PI)) * (2.*PI);
#endif
#if VSOP_GETRATE
/* convert millenium rate to day rate */
for (i = 3; i < 6; ++i) ret[i] /= VSOP_A1000;
#endif
#if VSOP_SPHERICAL
/* reduction from dynamical equinox of VSOP87 to FK5;
*/
if (prec < 5e-7) { /* 5e-7 rad = 0.1 arc seconds */
double L1, c1, s1;
L1 = ret[0] - degrad(13.97 * t[1] - 0.031 * t[2]);
c1 = cos(L1); s1 = sin(L1);
ret[0] += degrad(-0.09033 + 0.03916 * (c1 + s1) * tan(ret[1]))/3600.0;
ret[1] += degrad(0.03916 * (c1 - s1))/3600.0;
}
#endif
return (0);
}
Reference in New Issue View Git Blame Copy Permalink