1 #include "astro.h" 2 3 #define MAXE (.99) /* 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 50 ecc = elem.e; 51 if(ecc > MAXE) 52 ecc = MAXE; 53 incl = elem.i * radian; 54 node = (elem.o + 0.4593) * radian; 55 argp = (elem.w + elem.o + 0.4066) * radian; 56 mrad = elem.q / (1-ecc); 57 motion = .01720209895 * sqrt(1/(mrad*mrad*mrad))/radian; 58 anom = (eday - (elem.t - 2415020)) * motion * radian; 59 enom = anom + ecc*sin(anom); 60 61 do { 62 dele = (anom - enom + ecc * sin(enom)) / 63 (1 - ecc*cos(enom)); 64 enom += dele; 65 } while(fabs(dele) > 1.e-8); 66 67 vnom = 2*atan2( 68 sqrt((1+ecc)/(1-ecc))*sin(enom/2), 69 cos(enom/2)); 70 rad = mrad*(1-ecc*cos(enom)); 71 lambda = vnom + argp; 72 pturbl = 0; 73 lambda += pturbl*radsec; 74 pturbb = 0; 75 pturbr = 0; 76 77 /* 78 * reduce to the ecliptic 79 */ 80 nd = lambda - node; 81 lambda = node + atan2(sin(nd)*cos(incl),cos(nd)); 82 83 sl = sin(incl)*sin(nd) + pturbb*radsec; 84 beta = atan2(sl, sqrt(1-sl*sl)); 85 86 lograd = pturbr*2.30258509; 87 rad *= 1 + lograd; 88 89 motion *= radian*mrad*mrad/(rad*rad); 90 semi = 0; 91 92 mag = 5.47 + 6.1/2.303*log(rad); 93 94 helio(); 95 geo(); 96 } 97