sys/src/cmd/astro/comet.c

   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 }