1 files changed, 132 insertions, 0 deletions
diff --git a/src/cmd/astro/comet.c b/src/cmd/astro/comet.c
new file mode 100644
index 00000000..bfd11dc0
--- /dev/null
+++ b/src/cmd/astro/comet.c
@@ -0,0 +1,132 @@
+#include "astro.h"
+
+#define	MAXE	(.999)	/* cant do hyperbolas */
+
+void
+comet(void)
+{
+	double pturbl, pturbb, pturbr;
+	double lograd;
+	double dele, enom, vnom, nd, sl;
+
+	struct	elem
+	{
+		double	t;	/* time of perihelion */
+		double	q;	/* perihelion distance */
+		double	e;	/* eccentricity */
+		double	i;	/* inclination */
+		double	w;	/* argument of perihelion */
+		double	o;	/* longitude of ascending node */
+	} elem;
+
+/*	elem = (struct elem)
+	{
+		etdate(1990, 5, 19.293),
+		0.9362731,
+		0.6940149,
+		11.41096,
+		198.77059,
+		69.27130,
+	};	/* p/schwassmann-wachmann 3, 1989d */
+/*	elem = (struct elem)
+	{
+		etdate(1990, 4, 9.9761),
+		0.349957,
+		1.00038,
+		58.9596,
+		61.5546,
+		75.2132,
+	};	/* austin 3, 1989c */
+/*	elem = (struct elem)
+	{
+		etdate(1990, 10, 24.36),
+		0.9385,
+		1.00038,
+		131.62,
+		242.58,
+		138.57,
+	};	/* levy 6 , 1990c */
+/*	elem=(struct elem)
+	{
+		etdate(1996, 5, 1.3965),
+		0.230035,
+		0.999662,
+		124.9098,
+		130.2102,
+		188.0430,
+	};	/* C/1996 B2 (Hyakutake) */
+/*	elem=(struct elem)
+	{
+		etdate(1997, 4, 1.13413),
+		0.9141047,
+		0.9950989,
+		89.42932,
+		130.59066,
+		282.47069,
+	};	/*C/1995 O1 (Hale-Bopp) */
+/*	elem=(struct elem)
+	{
+		etdate(2000, 7, 26.1754),
+		0.765126,
+		0.999356,
+		149.3904,
+		151.0510,
+		83.1909,
+	};	/*C/1999 S4 (Linear) */
+	elem=(struct elem)
+	{
+		etdate(2002, 3, 18.9784),
+		0.5070601,
+		0.990111,
+		28.12106,
+		34.6666,
+		93.1206,
+	};	/*C/2002 C1 (Ikeya-Zhang) */
+
+	ecc = elem.e;
+	if(ecc > MAXE)
+		ecc = MAXE;
+	incl = elem.i * radian;
+	node = (elem.o + 0.4593) * radian;
+	argp = (elem.w + elem.o + 0.4066) * radian;
+	mrad = elem.q / (1-ecc);
+        motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian;
+	anom = (eday - (elem.t - 2415020)) * motion * 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, sqrt(1-sl*sl));
+
+	lograd = pturbr*2.30258509;
+	rad *= 1 + lograd;
+
+	motion *= radian*mrad*mrad/(rad*rad);
+	semi = 0;
+
+	mag = 5.47 + 6.1/2.303*log(rad);
+
+	helio();
+	geo();
+}