ref: 4308f6e7e6586e3a27aa46c4d54c434b8533b16d
dir: /sys/src/cmd/astro/comet.c/
#include "astro.h"
#define MAXE (.999) /* cant do hyperbolas */
void
comet(void)
{
double pturbl, pturbb, pturbr;
double lograd;
double dele, enom, vnom, nd, sl;
struct elem
{
double t; /* time of perihelion */
double q; /* perihelion distance */
double e; /* eccentricity */
double i; /* inclination */
double w; /* argument of perihelion */
double o; /* longitude of ascending node */
} elem;
/* elem = (struct elem)
{
etdate(1990, 5, 19.293),
0.9362731,
0.6940149,
11.41096,
198.77059,
69.27130,
}; /* p/schwassmann-wachmann 3, 1989d */
/* elem = (struct elem)
{
etdate(1990, 4, 9.9761),
0.349957,
1.00038,
58.9596,
61.5546,
75.2132,
}; /* austin 3, 1989c */
/* elem = (struct elem)
{
etdate(1990, 10, 24.36),
0.9385,
1.00038,
131.62,
242.58,
138.57,
}; /* levy 6 , 1990c */
/* elem=(struct elem)
{
etdate(1996, 5, 1.3965),
0.230035,
0.999662,
124.9098,
130.2102,
188.0430,
}; /* C/1996 B2 (Hyakutake) */
/* elem=(struct elem)
{
etdate(1997, 4, 1.13413),
0.9141047,
0.9950989,
89.42932,
130.59066,
282.47069,
}; /*C/1995 O1 (Hale-Bopp) */
/* elem=(struct elem)
{
etdate(2000, 7, 26.1754),
0.765126,
0.999356,
149.3904,
151.0510,
83.1909,
}; /*C/1999 S4 (Linear) */
elem=(struct elem)
{
etdate(2002, 3, 18.9784),
0.5070601,
0.990111,
28.12106,
34.6666,
93.1206,
}; /*C/2002 C1 (Ikeya-Zhang) */
ecc = elem.e;
if(ecc > MAXE)
ecc = MAXE;
incl = elem.i * radian;
node = (elem.o + 0.4593) * radian;
argp = (elem.w + elem.o + 0.4066) * radian;
mrad = elem.q / (1-ecc);
motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian;
anom = (eday - (elem.t - 2415020)) * motion * radian;
enom = anom + ecc*sin(anom);
do {
dele = (anom - enom + ecc * sin(enom)) /
(1 - ecc*cos(enom));
enom += dele;
} while(fabs(dele) > converge);
vnom = 2*atan2(
sqrt((1+ecc)/(1-ecc))*sin(enom/2),
cos(enom/2));
rad = mrad*(1-ecc*cos(enom));
lambda = vnom + argp;
pturbl = 0;
lambda += pturbl*radsec;
pturbb = 0;
pturbr = 0;
/*
* reduce to the ecliptic
*/
nd = lambda - node;
lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
sl = sin(incl)*sin(nd) + pturbb*radsec;
beta = atan2(sl, sqrt(1-sl*sl));
lograd = pturbr*2.30258509;
rad *= 1 + lograd;
motion *= radian*mrad*mrad/(rad*rad);
semi = 0;
mag = 5.47 + 6.1/2.303*log(rad);
helio();
geo();
}