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