mirror of
https://github.com/9fans/plan9port.git
synced 2025-01-15 11:20:03 +00:00
146 lines
2.5 KiB
C
146 lines
2.5 KiB
C
#include "astro.h"
|
|
|
|
#define MAXE (.999) /* cant do hyperbolas */
|
|
|
|
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 */
|
|
};
|
|
|
|
struct elem
|
|
mkelem(double t, double q, double e, double i, double w, double o)
|
|
{
|
|
struct elem el;
|
|
|
|
el.t = t;
|
|
el.q = q;
|
|
el.e = e;
|
|
el.i = i;
|
|
el.w = w;
|
|
el.o = o;
|
|
return el;
|
|
}
|
|
|
|
void
|
|
comet(void)
|
|
{
|
|
double pturbl, pturbb, pturbr;
|
|
double lograd;
|
|
double dele, enom, vnom, nd, sl;
|
|
struct elem 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 = mkelem(
|
|
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();
|
|
}
|