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