/*  * This file is part of the UCB release of Plan 9. It is subject to the license * terms in the LICENSE file found in the top-level directory of this * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No * part of the UCB release of Plan 9, including this file, may be copied, * modified, propagated, or distributed except according to the terms contained * in the LICENSE file. */ #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)}; voidplut(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();}