sys/src/cmd/astro/sat.c

   1 #include "astro.h"
   2
   3 void
   4 sat(void)
   5 {
   6         double pturbl, pturbb, pturbr;
   7         double lograd;
   8         double dele, enom, vnom, nd, sl;
   9
  10         double capj, capn, eye, comg, omg;
  11         double sb, su, cu, u, b, up;
  12         double sd, ca, sa;
  13
  14         ecc = .0558900 - .000347*capt;
  15         incl = 2.49256 - .0044*capt;
  16         node = 112.78364 + .87306*capt;
  17         argp = 91.08897 + 1.95917*capt;
  18         mrad = 9.538843;
  19         anom = 175.47630 + .03345972*eday - .56527*capt;
  20         motion = 120.4550/3600.;
  21
  22         incl *= radian;
  23         node *= radian;
  24         argp *= radian;
  25         anom = fmod(anom, 360.)*radian;
  26
  27         enom = anom + ecc*sin(anom);
  28         do {
  29                 dele = (anom - enom + ecc * sin(enom)) /
  30                         (1. - ecc*cos(enom));
  31                 enom += dele;
  32         } while(fabs(dele) > converge);
  33         vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  34                 cos(enom/2.));
  35         rad = mrad*(1. - ecc*cos(enom));
  36
  37         lambda = vnom + argp;
  38         pturbl = 0.;
  39         lambda += pturbl*radsec;
  40         pturbb = 0.;
  41         pturbr = 0.;
  42
  43 /*
  44  *      reduce to the ecliptic
  45  */
  46
  47         nd = lambda - node;
  48         lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
  49
  50         sl = sin(incl)*sin(nd) + pturbb*radsec;
  51         beta = atan2(sl, pyth(sl));
  52
  53         lograd = pturbr*2.30258509;
  54         rad *= 1. + lograd;
  55
  56
  57         lambda -= 1185.*radsec;
  58         beta -= 51.*radsec;
  59
  60         motion *= radian*mrad*mrad/(rad*rad);
  61         semi = 83.33;
  62
  63 /*
  64  *      here begins the computation of magnitude
  65  *      first find the geocentric equatorial coordinates of Saturn
  66  */
  67
  68         sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
  69                 sin(beta)*cos(obliq));
  70         sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
  71                 sin(beta)*sin(obliq));
  72         ca = rad*cos(beta)*cos(lambda);
  73         sd += zms;
  74         sa += yms;
  75         ca += xms;
  76         alpha = atan2(sa,ca);
  77         delta = atan2(sd,sqrt(sa*sa+ca*ca));
  78
  79 /*
  80  *      here are the necessary elements of Saturn's rings
  81  *      cf. Exp. Supp. p. 363ff.
  82  */
  83
  84         capj = 6.9056 - 0.4322*capt;
  85         capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
  86         eye = 28.0743 - 0.0128*capt;
  87         comg = 168.1179 + 1.3936*capt;
  88         omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
  89
  90         capj *= radian;
  91         capn *= radian;
  92         eye *= radian;
  93         comg *= radian;
  94         omg *= radian;
  95
  96 /*
  97  *      now find saturnicentric ring-plane coords of the earth
  98  */
  99
 100         sb = sin(capj)*cos(delta)*sin(alpha-capn) -
 101                 cos(capj)*sin(delta);
 102         su = cos(capj)*cos(delta)*sin(alpha-capn) +
 103                 sin(capj)*sin(delta);
 104         cu = cos(delta)*cos(alpha-capn);
 105         u = atan2(su,cu);
 106         b = atan2(sb,sqrt(su*su+cu*cu));
 107
 108 /*
 109  *      and then the saturnicentric ring-plane coords of the sun
 110  */
 111
 112         su = sin(eye)*sin(beta) +
 113                 cos(eye)*cos(beta)*sin(lambda-comg);
 114         cu = cos(beta)*cos(lambda-comg);
 115         up = atan2(su,cu);
 116
 117 /*
 118  *      at last, the magnitude
 119  */
 120
 121
 122         sb = sin(b);
 123         mag = -8.68 +2.52*fabs(up+omg-u)-
 124                 2.60*fabs(sb) + 1.25*(sb*sb);
 125
 126         helio();
 127         geo();
 128 }