13e12c5d1SDavid du Colombier #include "astro.h"
23e12c5d1SDavid du Colombier
33e12c5d1SDavid du Colombier void
sat(void)43e12c5d1SDavid du Colombier sat(void)
53e12c5d1SDavid du Colombier {
63e12c5d1SDavid du Colombier double pturbl, pturbb, pturbr;
73e12c5d1SDavid du Colombier double lograd;
83e12c5d1SDavid du Colombier double dele, enom, vnom, nd, sl;
93e12c5d1SDavid du Colombier
103e12c5d1SDavid du Colombier double capj, capn, eye, comg, omg;
113e12c5d1SDavid du Colombier double sb, su, cu, u, b, up;
123e12c5d1SDavid du Colombier double sd, ca, sa;
133e12c5d1SDavid du Colombier
143e12c5d1SDavid du Colombier ecc = .0558900 - .000347*capt;
153e12c5d1SDavid du Colombier incl = 2.49256 - .0044*capt;
163e12c5d1SDavid du Colombier node = 112.78364 + .87306*capt;
173e12c5d1SDavid du Colombier argp = 91.08897 + 1.95917*capt;
183e12c5d1SDavid du Colombier mrad = 9.538843;
193e12c5d1SDavid du Colombier anom = 175.47630 + .03345972*eday - .56527*capt;
203e12c5d1SDavid du Colombier motion = 120.4550/3600.;
213e12c5d1SDavid du Colombier
223e12c5d1SDavid du Colombier incl *= radian;
233e12c5d1SDavid du Colombier node *= radian;
243e12c5d1SDavid du Colombier argp *= radian;
253e12c5d1SDavid du Colombier anom = fmod(anom, 360.)*radian;
263e12c5d1SDavid du Colombier
273e12c5d1SDavid du Colombier enom = anom + ecc*sin(anom);
283e12c5d1SDavid du Colombier do {
293e12c5d1SDavid du Colombier dele = (anom - enom + ecc * sin(enom)) /
303e12c5d1SDavid du Colombier (1. - ecc*cos(enom));
313e12c5d1SDavid du Colombier enom += dele;
32*7dd7cddfSDavid du Colombier } while(fabs(dele) > converge);
333e12c5d1SDavid du Colombier vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
343e12c5d1SDavid du Colombier cos(enom/2.));
353e12c5d1SDavid du Colombier rad = mrad*(1. - ecc*cos(enom));
363e12c5d1SDavid du Colombier
373e12c5d1SDavid du Colombier lambda = vnom + argp;
383e12c5d1SDavid du Colombier pturbl = 0.;
393e12c5d1SDavid du Colombier lambda += pturbl*radsec;
403e12c5d1SDavid du Colombier pturbb = 0.;
413e12c5d1SDavid du Colombier pturbr = 0.;
423e12c5d1SDavid du Colombier
433e12c5d1SDavid du Colombier /*
443e12c5d1SDavid du Colombier * reduce to the ecliptic
453e12c5d1SDavid du Colombier */
463e12c5d1SDavid du Colombier
473e12c5d1SDavid du Colombier nd = lambda - node;
483e12c5d1SDavid du Colombier lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
493e12c5d1SDavid du Colombier
503e12c5d1SDavid du Colombier sl = sin(incl)*sin(nd) + pturbb*radsec;
513e12c5d1SDavid du Colombier beta = atan2(sl, pyth(sl));
523e12c5d1SDavid du Colombier
533e12c5d1SDavid du Colombier lograd = pturbr*2.30258509;
543e12c5d1SDavid du Colombier rad *= 1. + lograd;
553e12c5d1SDavid du Colombier
563e12c5d1SDavid du Colombier
573e12c5d1SDavid du Colombier lambda -= 1185.*radsec;
583e12c5d1SDavid du Colombier beta -= 51.*radsec;
593e12c5d1SDavid du Colombier
603e12c5d1SDavid du Colombier motion *= radian*mrad*mrad/(rad*rad);
613e12c5d1SDavid du Colombier semi = 83.33;
623e12c5d1SDavid du Colombier
633e12c5d1SDavid du Colombier /*
643e12c5d1SDavid du Colombier * here begins the computation of magnitude
653e12c5d1SDavid du Colombier * first find the geocentric equatorial coordinates of Saturn
663e12c5d1SDavid du Colombier */
673e12c5d1SDavid du Colombier
683e12c5d1SDavid du Colombier sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
693e12c5d1SDavid du Colombier sin(beta)*cos(obliq));
703e12c5d1SDavid du Colombier sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
713e12c5d1SDavid du Colombier sin(beta)*sin(obliq));
723e12c5d1SDavid du Colombier ca = rad*cos(beta)*cos(lambda);
733e12c5d1SDavid du Colombier sd += zms;
743e12c5d1SDavid du Colombier sa += yms;
753e12c5d1SDavid du Colombier ca += xms;
763e12c5d1SDavid du Colombier alpha = atan2(sa,ca);
773e12c5d1SDavid du Colombier delta = atan2(sd,sqrt(sa*sa+ca*ca));
783e12c5d1SDavid du Colombier
793e12c5d1SDavid du Colombier /*
803e12c5d1SDavid du Colombier * here are the necessary elements of Saturn's rings
813e12c5d1SDavid du Colombier * cf. Exp. Supp. p. 363ff.
823e12c5d1SDavid du Colombier */
833e12c5d1SDavid du Colombier
843e12c5d1SDavid du Colombier capj = 6.9056 - 0.4322*capt;
853e12c5d1SDavid du Colombier capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
863e12c5d1SDavid du Colombier eye = 28.0743 - 0.0128*capt;
873e12c5d1SDavid du Colombier comg = 168.1179 + 1.3936*capt;
883e12c5d1SDavid du Colombier omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
893e12c5d1SDavid du Colombier
903e12c5d1SDavid du Colombier capj *= radian;
913e12c5d1SDavid du Colombier capn *= radian;
923e12c5d1SDavid du Colombier eye *= radian;
933e12c5d1SDavid du Colombier comg *= radian;
943e12c5d1SDavid du Colombier omg *= radian;
953e12c5d1SDavid du Colombier
963e12c5d1SDavid du Colombier /*
973e12c5d1SDavid du Colombier * now find saturnicentric ring-plane coords of the earth
983e12c5d1SDavid du Colombier */
993e12c5d1SDavid du Colombier
1003e12c5d1SDavid du Colombier sb = sin(capj)*cos(delta)*sin(alpha-capn) -
1013e12c5d1SDavid du Colombier cos(capj)*sin(delta);
1023e12c5d1SDavid du Colombier su = cos(capj)*cos(delta)*sin(alpha-capn) +
1033e12c5d1SDavid du Colombier sin(capj)*sin(delta);
1043e12c5d1SDavid du Colombier cu = cos(delta)*cos(alpha-capn);
1053e12c5d1SDavid du Colombier u = atan2(su,cu);
1063e12c5d1SDavid du Colombier b = atan2(sb,sqrt(su*su+cu*cu));
1073e12c5d1SDavid du Colombier
1083e12c5d1SDavid du Colombier /*
1093e12c5d1SDavid du Colombier * and then the saturnicentric ring-plane coords of the sun
1103e12c5d1SDavid du Colombier */
1113e12c5d1SDavid du Colombier
1123e12c5d1SDavid du Colombier su = sin(eye)*sin(beta) +
1133e12c5d1SDavid du Colombier cos(eye)*cos(beta)*sin(lambda-comg);
1143e12c5d1SDavid du Colombier cu = cos(beta)*cos(lambda-comg);
1153e12c5d1SDavid du Colombier up = atan2(su,cu);
1163e12c5d1SDavid du Colombier
1173e12c5d1SDavid du Colombier /*
1183e12c5d1SDavid du Colombier * at last, the magnitude
1193e12c5d1SDavid du Colombier */
1203e12c5d1SDavid du Colombier
1213e12c5d1SDavid du Colombier
1223e12c5d1SDavid du Colombier sb = sin(b);
1233e12c5d1SDavid du Colombier mag = -8.68 +2.52*fabs(up+omg-u)-
1243e12c5d1SDavid du Colombier 2.60*fabs(sb) + 1.25*(sb*sb);
1253e12c5d1SDavid du Colombier
1263e12c5d1SDavid du Colombier helio();
1273e12c5d1SDavid du Colombier geo();
1283e12c5d1SDavid du Colombier }
129