#include "astro.h"

#define	MAXE	(.99)	/* 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 */

	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) > 1.e-8);

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