Subversion Repositories planix.SVN

#include "astro.h"

void
sat(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;

        ecc = .0558900 - .000347*capt;
        incl = 2.49256 - .0044*capt;
        node = 112.78364 + .87306*capt;
        argp = 91.08897 + 1.95917*capt;
        mrad = 9.538843;
        anom = 175.47630 + .03345972*eday - .56527*capt;
        motion = 120.4550/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();
}