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