1 cd5bae78 2004-04-21 devnull #include "astro.h"
  2 cd5bae78 2004-04-21 devnull 
  3 cd5bae78 2004-04-21 devnull #define	MAXE	(.999)	/* cant do hyperbolas */
  4 cd5bae78 2004-04-21 devnull 
  5 cd5bae78 2004-04-21 devnull void
  6 cd5bae78 2004-04-21 devnull comet(void)
  7 cd5bae78 2004-04-21 devnull {
  8 cd5bae78 2004-04-21 devnull 	double pturbl, pturbb, pturbr;
  9 cd5bae78 2004-04-21 devnull 	double lograd;
 10 cd5bae78 2004-04-21 devnull 	double dele, enom, vnom, nd, sl;
 11 cd5bae78 2004-04-21 devnull 
 12 cd5bae78 2004-04-21 devnull 	struct	elem
 13 cd5bae78 2004-04-21 devnull 	{
 14 cd5bae78 2004-04-21 devnull 		double	t;	/* time of perihelion */
 15 cd5bae78 2004-04-21 devnull 		double	q;	/* perihelion distance */
 16 cd5bae78 2004-04-21 devnull 		double	e;	/* eccentricity */
 17 cd5bae78 2004-04-21 devnull 		double	i;	/* inclination */
 18 cd5bae78 2004-04-21 devnull 		double	w;	/* argument of perihelion */
 19 cd5bae78 2004-04-21 devnull 		double	o;	/* longitude of ascending node */
 20 cd5bae78 2004-04-21 devnull 	} elem;
 21 cd5bae78 2004-04-21 devnull 
 22 cd5bae78 2004-04-21 devnull /*	elem = (struct elem)
 23 cd5bae78 2004-04-21 devnull 	{
 24 cd5bae78 2004-04-21 devnull 		etdate(1990, 5, 19.293),
 25 cd5bae78 2004-04-21 devnull 		0.9362731,
 26 cd5bae78 2004-04-21 devnull 		0.6940149,
 27 cd5bae78 2004-04-21 devnull 		11.41096,
 28 cd5bae78 2004-04-21 devnull 		198.77059,
 29 cd5bae78 2004-04-21 devnull 		69.27130,
 30 cd5bae78 2004-04-21 devnull 	};	/* p/schwassmann-wachmann 3, 1989d */
 31 cd5bae78 2004-04-21 devnull /*	elem = (struct elem)
 32 cd5bae78 2004-04-21 devnull 	{
 33 cd5bae78 2004-04-21 devnull 		etdate(1990, 4, 9.9761),
 34 cd5bae78 2004-04-21 devnull 		0.349957,
 35 cd5bae78 2004-04-21 devnull 		1.00038,
 36 cd5bae78 2004-04-21 devnull 		58.9596,
 37 cd5bae78 2004-04-21 devnull 		61.5546,
 38 cd5bae78 2004-04-21 devnull 		75.2132,
 39 cd5bae78 2004-04-21 devnull 	};	/* austin 3, 1989c */
 40 cd5bae78 2004-04-21 devnull /*	elem = (struct elem)
 41 cd5bae78 2004-04-21 devnull 	{
 42 cd5bae78 2004-04-21 devnull 		etdate(1990, 10, 24.36),
 43 cd5bae78 2004-04-21 devnull 		0.9385,
 44 cd5bae78 2004-04-21 devnull 		1.00038,
 45 cd5bae78 2004-04-21 devnull 		131.62,
 46 cd5bae78 2004-04-21 devnull 		242.58,
 47 cd5bae78 2004-04-21 devnull 		138.57,
 48 cd5bae78 2004-04-21 devnull 	};	/* levy 6 , 1990c */
 49 cd5bae78 2004-04-21 devnull /*	elem=(struct elem)
 50 cd5bae78 2004-04-21 devnull 	{
 51 cd5bae78 2004-04-21 devnull 		etdate(1996, 5, 1.3965),
 52 cd5bae78 2004-04-21 devnull 		0.230035,
 53 cd5bae78 2004-04-21 devnull 		0.999662,
 54 cd5bae78 2004-04-21 devnull 		124.9098,
 55 cd5bae78 2004-04-21 devnull 		130.2102,
 56 cd5bae78 2004-04-21 devnull 		188.0430,
 57 cd5bae78 2004-04-21 devnull 	};	/* C/1996 B2 (Hyakutake) */
 58 cd5bae78 2004-04-21 devnull /*	elem=(struct elem)
 59 cd5bae78 2004-04-21 devnull 	{
 60 cd5bae78 2004-04-21 devnull 		etdate(1997, 4, 1.13413),
 61 cd5bae78 2004-04-21 devnull 		0.9141047,
 62 cd5bae78 2004-04-21 devnull 		0.9950989,
 63 cd5bae78 2004-04-21 devnull 		89.42932,
 64 cd5bae78 2004-04-21 devnull 		130.59066,
 65 cd5bae78 2004-04-21 devnull 		282.47069,
 66 cd5bae78 2004-04-21 devnull 	};	/*C/1995 O1 (Hale-Bopp) */
 67 cd5bae78 2004-04-21 devnull /*	elem=(struct elem)
 68 cd5bae78 2004-04-21 devnull 	{
 69 cd5bae78 2004-04-21 devnull 		etdate(2000, 7, 26.1754),
 70 cd5bae78 2004-04-21 devnull 		0.765126,
 71 cd5bae78 2004-04-21 devnull 		0.999356,
 72 cd5bae78 2004-04-21 devnull 		149.3904,
 73 cd5bae78 2004-04-21 devnull 		151.0510,
 74 cd5bae78 2004-04-21 devnull 		83.1909,
 75 cd5bae78 2004-04-21 devnull 	};	/*C/1999 S4 (Linear) */
 76 cd5bae78 2004-04-21 devnull 	elem=(struct elem)
 77 cd5bae78 2004-04-21 devnull 	{
 78 cd5bae78 2004-04-21 devnull 		etdate(2002, 3, 18.9784),
 79 cd5bae78 2004-04-21 devnull 		0.5070601,
 80 cd5bae78 2004-04-21 devnull 		0.990111,
 81 cd5bae78 2004-04-21 devnull 		28.12106,
 82 cd5bae78 2004-04-21 devnull 		34.6666,
 83 cd5bae78 2004-04-21 devnull 		93.1206,
 84 cd5bae78 2004-04-21 devnull 	};	/*C/2002 C1 (Ikeya-Zhang) */
 85 cd5bae78 2004-04-21 devnull 
 86 cd5bae78 2004-04-21 devnull 	ecc = elem.e;
 87 cd5bae78 2004-04-21 devnull 	if(ecc > MAXE)
 88 cd5bae78 2004-04-21 devnull 		ecc = MAXE;
 89 cd5bae78 2004-04-21 devnull 	incl = elem.i * radian;
 90 cd5bae78 2004-04-21 devnull 	node = (elem.o + 0.4593) * radian;
 91 cd5bae78 2004-04-21 devnull 	argp = (elem.w + elem.o + 0.4066) * radian;
 92 cd5bae78 2004-04-21 devnull 	mrad = elem.q / (1-ecc);
 93 cd5bae78 2004-04-21 devnull         motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian;
 94 cd5bae78 2004-04-21 devnull 	anom = (eday - (elem.t - 2415020)) * motion * radian;
 95 cd5bae78 2004-04-21 devnull 	enom = anom + ecc*sin(anom);
 96 cd5bae78 2004-04-21 devnull 
 97 cd5bae78 2004-04-21 devnull 	do {
 98 cd5bae78 2004-04-21 devnull 		dele = (anom - enom + ecc * sin(enom)) /
 99 cd5bae78 2004-04-21 devnull 			(1 - ecc*cos(enom));
100 cd5bae78 2004-04-21 devnull 		enom += dele;
101 cd5bae78 2004-04-21 devnull 	} while(fabs(dele) > converge);
102 cd5bae78 2004-04-21 devnull 
103 cd5bae78 2004-04-21 devnull 	vnom = 2*atan2(
104 cd5bae78 2004-04-21 devnull 		sqrt((1+ecc)/(1-ecc))*sin(enom/2),
105 cd5bae78 2004-04-21 devnull 		cos(enom/2));
106 cd5bae78 2004-04-21 devnull 	rad = mrad*(1-ecc*cos(enom));
107 cd5bae78 2004-04-21 devnull 	lambda = vnom + argp;
108 cd5bae78 2004-04-21 devnull 	pturbl = 0;
109 cd5bae78 2004-04-21 devnull 	lambda += pturbl*radsec;
110 cd5bae78 2004-04-21 devnull 	pturbb = 0;
111 cd5bae78 2004-04-21 devnull 	pturbr = 0;
112 cd5bae78 2004-04-21 devnull 
113 cd5bae78 2004-04-21 devnull /*
114 cd5bae78 2004-04-21 devnull  *	reduce to the ecliptic
115 cd5bae78 2004-04-21 devnull  */
116 cd5bae78 2004-04-21 devnull 	nd = lambda - node;
117 cd5bae78 2004-04-21 devnull 	lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
118 cd5bae78 2004-04-21 devnull 
119 cd5bae78 2004-04-21 devnull 	sl = sin(incl)*sin(nd) + pturbb*radsec;
120 cd5bae78 2004-04-21 devnull 	beta = atan2(sl, sqrt(1-sl*sl));
121 cd5bae78 2004-04-21 devnull 
122 cd5bae78 2004-04-21 devnull 	lograd = pturbr*2.30258509;
123 cd5bae78 2004-04-21 devnull 	rad *= 1 + lograd;
124 cd5bae78 2004-04-21 devnull 
125 cd5bae78 2004-04-21 devnull 	motion *= radian*mrad*mrad/(rad*rad);
126 cd5bae78 2004-04-21 devnull 	semi = 0;
127 cd5bae78 2004-04-21 devnull 
128 cd5bae78 2004-04-21 devnull 	mag = 5.47 + 6.1/2.303*log(rad);
129 cd5bae78 2004-04-21 devnull 
130 cd5bae78 2004-04-21 devnull 	helio();
131 cd5bae78 2004-04-21 devnull 	geo();
132 cd5bae78 2004-04-21 devnull }