3 #define MAXE (.999) /* cant do hyperbolas */
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 */
16 mkelem(double t, double q, double e, double i, double w, double o)
32 double pturbl, pturbb, pturbr;
34 double dele, enom, vnom, nd, sl;
37 /* elem = (struct elem)
39 etdate(1990, 5, 19.293),
45 }; /* p/schwassmann-wachmann 3, 1989d */
46 /* elem = (struct elem)
48 etdate(1990, 4, 9.9761),
54 }; /* austin 3, 1989c */
55 /* elem = (struct elem)
57 etdate(1990, 10, 24.36),
63 }; /* levy 6 , 1990c */
66 etdate(1996, 5, 1.3965),
72 }; /* C/1996 B2 (Hyakutake) */
75 etdate(1997, 4, 1.13413),
81 }; /*C/1995 O1 (Hale-Bopp) */
84 etdate(2000, 7, 26.1754),
90 }; /*C/1999 S4 (Linear) */
92 etdate(2002, 3, 18.9784),
98 ); /*C/2002 C1 (Ikeya-Zhang) */
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);
112 dele = (anom - enom + ecc * sin(enom)) /
115 } while(fabs(dele) > converge);
118 sqrt((1+ecc)/(1-ecc))*sin(enom/2),
120 rad = mrad*(1-ecc*cos(enom));
121 lambda = vnom + argp;
123 lambda += pturbl*radsec;
128 * reduce to the ecliptic
131 lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
133 sl = sin(incl)*sin(nd) + pturbb*radsec;
134 beta = atan2(sl, sqrt(1-sl*sl));
136 lograd = pturbr*2.30258509;
139 motion *= radian*mrad*mrad/(rad*rad);
142 mag = 5.47 + 6.1/2.303*log(rad);