1 #include "astro.h"
2
3 #define MAXE (.999) /* cant do hyperbolas */
4
5 void
comet(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 }
133