ref: d2b3e2f55a52df45c8e9c0b043763809f0001e69
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(); }