#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();
}