#include "astro.h"

#define	MAXE	(.999)	/* cant do hyperbolas */

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 */
};

struct elem
mkelem(double t, double q, double e, double i, double w, double o)
{
	struct elem el;

	el.t = t;
	el.q = q;
	el.e = e;
	el.i = i;
	el.w = w;
	el.o = o;
	return el;
}

void
comet(void)
{
	double pturbl, pturbb, pturbr;
	double lograd;
	double dele, enom, vnom, nd, sl;
	struct elem 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 = mkelem(
		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();
}