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

}