mirror of https://github.com/XEphem/XEphem.git
228 lines
7.2 KiB
C
228 lines
7.2 KiB
C
/* rewritten for Bureau des Longitude theories by Bretagnon and Chapront
|
|
* Michael Sternberg <sternberg@physik.tu-chemnitz.de>
|
|
*/
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
#include "astro.h"
|
|
#include "vsop87.h"
|
|
#include "chap95.h"
|
|
|
|
static void pluto_ell (double mj, double *ret);
|
|
static void chap_trans (double mj, double *ret);
|
|
static void planpos (double mj, int obj, double prec, double *ret);
|
|
|
|
/* coordinate transformation
|
|
* from:
|
|
* J2000.0 rectangular equatoreal ret[{0,1,2}] = {x,y,z}
|
|
* to:
|
|
* mean equinox of date spherical ecliptical ret[{0,1,2}] = {l,b,r}
|
|
*/
|
|
static void
|
|
chap_trans (
|
|
double mj, /* destination epoch */
|
|
double *ret) /* vector to be transformed _IN PLACE_ */
|
|
{
|
|
double ra, dec, r, eps;
|
|
double sr, cr, sd, cd, se, ce;
|
|
|
|
cartsph(ret[0], ret[1], ret[2], &ra, &dec, &r);
|
|
precess(J2000, mj, &ra, &dec);
|
|
obliquity(mj, &eps);
|
|
sr = sin(ra); cr = cos(ra);
|
|
sd = sin(dec); cd = cos(dec);
|
|
se = sin(eps); ce = cos(eps);
|
|
ret[0] = atan2( sr * ce + sd/cd * se, cr); /* long */
|
|
ret[1] = asin( sd * ce - cd * se * sr); /* lat */
|
|
ret[2] = r; /* radius */
|
|
}
|
|
|
|
/* low precision ecliptic coordinates of Pluto from mean orbit.
|
|
* Only for sake of completeness outside available perturbation theories.
|
|
*/
|
|
static void
|
|
pluto_ell (
|
|
double mj, /* epoch */
|
|
double *ret) /* ecliptic coordinates {l,b,r} at equinox of date */
|
|
{
|
|
/* Mean orbital elements of Pluto.
|
|
* Astronomical Almanac 2020 page G3.
|
|
*/
|
|
double a = 39.789, /* semimajor axis, au */
|
|
e = 0.252, /* excentricity */
|
|
inc0 = 17.097, /* inclination, deg */
|
|
Om0 = 110.297, /* long asc node, deg */
|
|
omeg0 = 115.058, /* arg of perihel, deg */
|
|
mjp = 2459000.5 - MJD0, /* epoch of perihel */
|
|
mjeq = J2000, /* equinox of elements */
|
|
n = 0.0039; /* daily motion, deg */
|
|
|
|
double inc, Om, omeg; /* orbital elements at epoch of date */
|
|
double ma, ea, nu; /* mean, excentric and true anomaly */
|
|
double lo, slo, clo; /* longitude in orbit from asc node */
|
|
|
|
reduce_elements(mjeq, mj, degrad(inc0), degrad(omeg0), degrad(Om0),
|
|
&inc, &omeg, &Om);
|
|
ma = degrad((mj - mjp) * n);
|
|
anomaly(ma, e, &nu, &ea);
|
|
ret[2] = a * (1.0 - e*cos(ea)); /* r */
|
|
lo = omeg + nu;
|
|
slo = sin(lo);
|
|
clo = cos(lo);
|
|
ret[1] = asin(slo * sin(inc)); /* b */
|
|
ret[0] = atan2(slo * cos(inc), clo) + Om; /* l */
|
|
}
|
|
|
|
/*************************************************************/
|
|
|
|
/* geometric heliocentric position of planet, mean ecliptic of date
|
|
* (not corrected for light-time)
|
|
*/
|
|
static void
|
|
planpos (double mj, int obj, double prec, double *ret)
|
|
{
|
|
if (mj >= CHAP_BEGIN && mj <= CHAP_END) {
|
|
if (obj >= JUPITER) { /* prefer Chapront */
|
|
chap95(mj, obj, prec, ret);
|
|
chap_trans (mj, ret);
|
|
} else { /* VSOP for inner planets */
|
|
vsop87(mj, obj, prec, ret);
|
|
}
|
|
} else { /* outside Chapront time: */
|
|
if (obj != PLUTO) { /* VSOP for all but Pluto */
|
|
vsop87(mj, obj, prec, ret);
|
|
} else { /* Pluto mean elliptic orbit */
|
|
pluto_ell(mj, ret);
|
|
}
|
|
}
|
|
}
|
|
|
|
/*************************************************************/
|
|
|
|
/* visual elements of planets
|
|
* [planet][0] = angular size at 1 AU
|
|
* [planet][1] = magnitude at 1 AU from sun and earth and 0 deg phase angle
|
|
* [planet][2] = A
|
|
* [planet][3] = B
|
|
* [planet][4] = C
|
|
* where mag correction = A*(i/100) + B*(i/100)^2 + C*(i/100)^3
|
|
* i = angle between sun and earth from planet, degrees
|
|
* from Explanatory Supplement, 1992
|
|
*/
|
|
static double vis_elements[8][5] = {
|
|
/* Mercury */ { 6.74, -0.36, 3.8, -2.73, 2.00},
|
|
/* Venus */ { 16.92, -4.29, 0.09, 2.39, -.65},
|
|
/* Mars */ { 9.36, -1.52, 1.60, 0., 0.},
|
|
/* Jupiter */ { 196.74, -9.25, 0.50, 0., 0.},
|
|
/* Saturn */ { 165.6, -8.88, 4.40, 0., 0.},
|
|
/* Uranus */ { 65.8, -7.19, 0.28, 0., 0.},
|
|
/* Neptune */ { 62.2, -6.87, 0., 0., 0.},
|
|
/* Pluto */ { 8.2, -1.01, 4.1, 0., 0.}
|
|
};
|
|
|
|
/* given a modified Julian date, mj, and a planet, p, find:
|
|
* lpd0: heliocentric longitude,
|
|
* psi0: heliocentric latitude,
|
|
* rp0: distance from the sun to the planet,
|
|
* rho0: distance from the Earth to the planet,
|
|
* none corrected for light time, ie, they are the true values for the
|
|
* given instant.
|
|
* lam: geocentric ecliptic longitude,
|
|
* bet: geocentric ecliptic latitude,
|
|
* each corrected for light time, ie, they are the apparent values as
|
|
* seen from the center of the Earth for the given instant.
|
|
* dia: angular diameter in arcsec at 1 AU,
|
|
* mag: visual magnitude
|
|
*
|
|
* all angles are in radians, all distances in AU.
|
|
*
|
|
* corrections for nutation and abberation must be made by the caller. The RA
|
|
* and DEC calculated from the fully-corrected ecliptic coordinates are then
|
|
* the apparent geocentric coordinates. Further corrections can be made, if
|
|
* required, for atmospheric refraction and geocentric parallax.
|
|
*/
|
|
void
|
|
plans (double mj, PLCode p, double *lpd0, double *psi0, double *rp0,
|
|
double *rho0, double *lam, double *bet, double *dia, double *mag)
|
|
{
|
|
static double lastmj = -10000;
|
|
static double lsn, bsn, rsn; /* geocentric coords of sun */
|
|
static double xsn, ysn, zsn; /* cartesian " */
|
|
double lp, bp, rp; /* heliocentric coords of planet */
|
|
double xp, yp, zp, rho; /* rect. coords and geocentric dist. */
|
|
double dt; /* light time */
|
|
double *vp; /* vis_elements[p] */
|
|
double ci, i; /* sun/earth angle: cos, degrees */
|
|
int pass;
|
|
|
|
/* get sun cartesian; needed only once at mj */
|
|
if (mj != lastmj) {
|
|
sunpos (mj, &lsn, &rsn, &bsn);
|
|
sphcart (lsn, bsn, rsn, &xsn, &ysn, &zsn);
|
|
lastmj = mj;
|
|
}
|
|
|
|
/* first find the true position of the planet at mj.
|
|
* then repeat a second time for a slightly different time based
|
|
* on the position found in the first pass to account for light-travel
|
|
* time.
|
|
*/
|
|
dt = 0.0;
|
|
for (pass = 0; pass < 2; pass++) {
|
|
double ret[6];
|
|
|
|
/* get spherical coordinates of planet from precision routines,
|
|
* retarded for light time in second pass;
|
|
* alternative option: vsop allows calculating rates.
|
|
*/
|
|
planpos(mj - dt, p, 0.0, ret);
|
|
|
|
lp = ret[0];
|
|
bp = ret[1];
|
|
rp = ret[2];
|
|
|
|
sphcart (lp, bp, rp, &xp, &yp, &zp);
|
|
cartsph (xp + xsn, yp + ysn, zp + zsn, lam, bet, &rho);
|
|
|
|
if (pass == 0) {
|
|
/* save heliocentric coordinates at first pass since, being
|
|
* true, they are NOT to be corrected for light-travel time.
|
|
*/
|
|
*lpd0 = lp;
|
|
range (lpd0, 2.*PI);
|
|
*psi0 = bp;
|
|
*rp0 = rp;
|
|
*rho0 = rho;
|
|
}
|
|
|
|
/* when we view a planet we see it in the position it occupied
|
|
* dt days ago, where rho is the distance between it and earth,
|
|
* in AU. use this as the new time for the next pass.
|
|
*/
|
|
dt = rho * 5.7755183e-3;
|
|
}
|
|
|
|
vp = vis_elements[p];
|
|
*dia = vp[0];
|
|
|
|
/* solve plane triangle, assume sun/earth dist == 1 */
|
|
ci = (rp*rp + rho*rho - 1)/(2*rp*rho);
|
|
|
|
/* expl supp equation for mag */
|
|
if (ci < -1) ci = -1;
|
|
if (ci > 1) ci = 1;
|
|
i = raddeg(acos(ci))/100.;
|
|
*mag = vp[1] + 5*log10(rho*rp) + i*(vp[2] + i*(vp[3] + i*vp[4]));
|
|
|
|
/* rings contribution if SATURN */
|
|
if (p == SATURN) {
|
|
double et, st, set;
|
|
satrings (bp, lp, rp, lsn+PI, rsn, mj+MJD0, &et, &st);
|
|
set = sin(fabs(et));
|
|
*mag += (-2.60 + 1.25*set)*set;
|
|
}
|
|
}
|
|
|
|
/* For RCS Only -- Do Not Edit */
|
|
static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: plans.c,v $ $Date: 2003/11/15 04:07:36 $ $Revision: 1.5 $ $Name: $"};
|