1 #include "astro.h" 2 3 #define MAXE (.999) /* cant do hyperbolas */ 4 5 void 6 comet(void) 7 { 8 double pturbl, pturbb, pturbr; 9 double lograd; 10 double dele, enom, vnom, nd, sl; 11 12 struct elem 13 { 14 double t; /* time of perihelion */ 15 double q; /* perihelion distance */ 16 double e; /* eccentricity */ 17 double i; /* inclination */ 18 double w; /* argument of perihelion */ 19 double o; /* longitude of ascending node */ 20 } elem; 21 22 /* elem = (struct elem) 23 { 24 etdate(1990, 5, 19.293), 25 0.9362731, 26 0.6940149, 27 11.41096, 28 198.77059, 29 69.27130, 30 }; /* p/schwassmann-wachmann 3, 1989d */ 31 /* elem = (struct elem) 32 { 33 etdate(1990, 4, 9.9761), 34 0.349957, 35 1.00038, 36 58.9596, 37 61.5546, 38 75.2132, 39 }; /* austin 3, 1989c */ 40 /* elem = (struct elem) 41 { 42 etdate(1990, 10, 24.36), 43 0.9385, 44 1.00038, 45 131.62, 46 242.58, 47 138.57, 48 }; /* levy 6 , 1990c */ 49 /* elem=(struct elem) 50 { 51 etdate(1996, 5, 1.3965), 52 0.230035, 53 0.999662, 54 124.9098, 55 130.2102, 56 188.0430, 57 }; /* C/1996 B2 (Hyakutake) */ 58 /* elem=(struct elem) 59 { 60 etdate(1997, 4, 1.13413), 61 0.9141047, 62 0.9950989, 63 89.42932, 64 130.59066, 65 282.47069, 66 }; /*C/1995 O1 (Hale-Bopp) */ 67 /* elem=(struct elem) 68 { 69 etdate(2000, 7, 26.1754), 70 0.765126, 71 0.999356, 72 149.3904, 73 151.0510, 74 83.1909, 75 }; /*C/1999 S4 (Linear) */ 76 elem=(struct elem) 77 { 78 etdate(2002, 3, 18.9784), 79 0.5070601, 80 0.990111, 81 28.12106, 82 34.6666, 83 93.1206, 84 }; /*C/2002 C1 (Ikeya-Zhang) */ 85 86 ecc = elem.e; 87 if(ecc > MAXE) 88 ecc = MAXE; 89 incl = elem.i * radian; 90 node = (elem.o + 0.4593) * radian; 91 argp = (elem.w + elem.o + 0.4066) * radian; 92 mrad = elem.q / (1-ecc); 93 motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian; 94 anom = (eday - (elem.t - 2415020)) * motion * radian; 95 enom = anom + ecc*sin(anom); 96 97 do { 98 dele = (anom - enom + ecc * sin(enom)) / 99 (1 - ecc*cos(enom)); 100 enom += dele; 101 } while(fabs(dele) > converge); 102 103 vnom = 2*atan2( 104 sqrt((1+ecc)/(1-ecc))*sin(enom/2), 105 cos(enom/2)); 106 rad = mrad*(1-ecc*cos(enom)); 107 lambda = vnom + argp; 108 pturbl = 0; 109 lambda += pturbl*radsec; 110 pturbb = 0; 111 pturbr = 0; 112 113 /* 114 * reduce to the ecliptic 115 */ 116 nd = lambda - node; 117 lambda = node + atan2(sin(nd)*cos(incl),cos(nd)); 118 119 sl = sin(incl)*sin(nd) + pturbb*radsec; 120 beta = atan2(sl, sqrt(1-sl*sl)); 121 122 lograd = pturbr*2.30258509; 123 rad *= 1 + lograd; 124 125 motion *= radian*mrad*mrad/(rad*rad); 126 semi = 0; 127 128 mag = 5.47 + 6.1/2.303*log(rad); 129 130 helio(); 131 geo(); 132 } 133