123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132 |
- #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();
- }
|