3 #define MAXE (.999) /* cant do hyperbolas */
8 double pturbl, pturbb, pturbr;
10 double dele, enom, vnom, nd, sl;
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 */
22 /* elem = (struct elem)
24 etdate(1990, 5, 19.293),
30 }; /* p/schwassmann-wachmann 3, 1989d */
31 /* elem = (struct elem)
33 etdate(1990, 4, 9.9761),
39 }; /* austin 3, 1989c */
40 /* elem = (struct elem)
42 etdate(1990, 10, 24.36),
48 }; /* levy 6 , 1990c */
51 etdate(1996, 5, 1.3965),
57 }; /* C/1996 B2 (Hyakutake) */
60 etdate(1997, 4, 1.13413),
66 }; /*C/1995 O1 (Hale-Bopp) */
69 etdate(2000, 7, 26.1754),
75 }; /*C/1999 S4 (Linear) */
78 etdate(2002, 3, 18.9784),
84 }; /*C/2002 C1 (Ikeya-Zhang) */
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);
98 dele = (anom - enom + ecc * sin(enom)) /
101 } while(fabs(dele) > converge);
104 sqrt((1+ecc)/(1-ecc))*sin(enom/2),
106 rad = mrad*(1-ecc*cos(enom));
107 lambda = vnom + argp;
109 lambda += pturbl*radsec;
114 * reduce to the ecliptic
117 lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
119 sl = sin(incl)*sin(nd) + pturbb*radsec;
120 beta = atan2(sl, sqrt(1-sl*sl));
122 lograd = pturbr*2.30258509;
125 motion *= radian*mrad*mrad/(rad*rad);
128 mag = 5.47 + 6.1/2.303*log(rad);