blob: 6c31ed5894721ee971448b24cf60821bf483a313 [file] [log] [blame]
#include "astro.h"
static double elem[] =
{
36525.0, /* [0] eday of epoc */
39.48168677, /* [1] semi major axis (au) */
0.24880766, /* [2] eccentricity */
17.14175, /* [3] inclination (deg) */
110.30347, /* [4] longitude of the ascending node (deg) */
224.06676, /* [5] longitude of perihelion (deg) */
238.92881, /* [6] mean longitude (deg) */
-0.00076912, /* [1+6] (au/julian century) */
0.00006465, /* [2+6] (e/julian century) */
11.07, /* [3+6] (arcsec/julian century) */
-37.33, /* [4+6] (arcsec/julian century) */
-132.25, /* [5+6] (arcsec/julian century) */
522747.90, /* [6+6] (arcsec/julian century) */
};
void
plut(void)
{
double pturbl, pturbb, pturbr;
double lograd;
double dele, enom, vnom, nd, sl;
double capj, capn, eye, comg, omg;
double sb, su, cu, u, b, up;
double sd, ca, sa;
double cy;
cy = (eday - elem[0]) / 36525.; /* per julian century */
mrad = elem[1] + elem[1+6]*cy;
ecc = elem[2] + elem[2+6]*cy;
cy = cy / 3600; /* arcsec/deg per julian century */
incl = elem[3] + elem[3+6]*cy;
node = elem[4] + elem[4+6]*cy;
argp = elem[5] + elem[5+6]*cy;
anom = elem[6] + elem[6+6]*cy - argp;
motion = elem[6+6] / 36525. / 3600;
incl *= radian;
node *= radian;
argp *= radian;
anom = fmod(anom,360.)*radian;
enom = anom + ecc*sin(anom);
do {
dele = (anom - enom + ecc * sin(enom)) /
(1. - ecc*cos(enom));
enom += dele;
} while(fabs(dele) > converge);
vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
cos(enom/2.));
rad = mrad*(1. - ecc*cos(enom));
lambda = vnom + argp;
pturbl = 0.;
lambda += pturbl*radsec;
pturbb = 0.;
pturbr = 0.;
/*
* reduce to the ecliptic
*/
nd = lambda - node;
lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
sl = sin(incl)*sin(nd) + pturbb*radsec;
beta = atan2(sl, pyth(sl));
lograd = pturbr*2.30258509;
rad *= 1. + lograd;
lambda -= 1185.*radsec;
beta -= 51.*radsec;
motion *= radian*mrad*mrad/(rad*rad);
semi = 83.33;
/*
* here begins the computation of magnitude
* first find the geocentric equatorial coordinates of Saturn
*/
sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
sin(beta)*cos(obliq));
sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
sin(beta)*sin(obliq));
ca = rad*cos(beta)*cos(lambda);
sd += zms;
sa += yms;
ca += xms;
alpha = atan2(sa,ca);
delta = atan2(sd,sqrt(sa*sa+ca*ca));
/*
* here are the necessary elements of Saturn's rings
* cf. Exp. Supp. p. 363ff.
*/
capj = 6.9056 - 0.4322*capt;
capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
eye = 28.0743 - 0.0128*capt;
comg = 168.1179 + 1.3936*capt;
omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
capj *= radian;
capn *= radian;
eye *= radian;
comg *= radian;
omg *= radian;
/*
* now find saturnicentric ring-plane coords of the earth
*/
sb = sin(capj)*cos(delta)*sin(alpha-capn) -
cos(capj)*sin(delta);
su = cos(capj)*cos(delta)*sin(alpha-capn) +
sin(capj)*sin(delta);
cu = cos(delta)*cos(alpha-capn);
u = atan2(su,cu);
b = atan2(sb,sqrt(su*su+cu*cu));
/*
* and then the saturnicentric ring-plane coords of the sun
*/
su = sin(eye)*sin(beta) +
cos(eye)*cos(beta)*sin(lambda-comg);
cu = cos(beta)*cos(lambda-comg);
up = atan2(su,cu);
/*
* at last, the magnitude
*/
sb = sin(b);
mag = -8.68 +2.52*fabs(up+omg-u)-
2.60*fabs(sb) + 1.25*(sb*sb);
helio();
geo();
}