Astro with some minor changes to placate Unix.

1 files changed, 154 insertions, 0 deletions

diff --git a/src/cmd/astro/plut.c b/src/cmd/astro/plut.c
new file mode 100644
index 00000000..36702a80
--- /dev/null
+++ b/src/cmd/astro/plut.c
@@ -0,0 +1,154 @@
+#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();
+}