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