/*  * 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" #define	MAXE	(.999)	/* cant do hyperbolas */ voidcomet(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();}