sys/src/cmd/astro/nept.c

   1 #include "astro.h"
   2
   3 static  double  elem[] =
   4 {
   5         36525.0,                // [0] eday of epoc
   6
   7         30.06896348,            // [1] semi major axis (au)
   8         0.00858587,             // [2] eccentricity
   9         1.76917,                // [3] inclination (deg)
  10         131.72169,              // [4] longitude of the ascending node (deg)
  11         44.97135,               // [5] longitude of perihelion (deg)
  12         304.88003,              // [6] mean longitude (deg)
  13
  14         -0.00125196,            // [1+6] (au/julian century)
  15         0.0000251,              // [2+6] (e/julian century)
  16         -3.64,                  // [3+6] (arcsec/julian century)
  17         -151.25,                // [4+6] (arcsec/julian century)
  18         -844.43,                // [5+6] (arcsec/julian century)
  19         786449.21,              // [6+6] (arcsec/julian century)
  20 };
  21
  22 void
  23 nept(void)
  24 {
  25         double pturbl, pturbb, pturbr;
  26         double lograd;
  27         double dele, enom, vnom, nd, sl;
  28
  29         double capj, capn, eye, comg, omg;
  30         double sb, su, cu, u, b, up;
  31         double sd, ca, sa;
  32
  33         double cy;
  34
  35         cy = (eday - elem[0]) / 36525.;         // per julian century
  36
  37         mrad = elem[1] + elem[1+6]*cy;
  38         ecc = elem[2] + elem[2+6]*cy;
  39
  40         cy = cy / 3600;                         // arcsec/deg per julian century
  41         incl = elem[3] + elem[3+6]*cy;
  42         node = elem[4] + elem[4+6]*cy;
  43         argp = elem[5] + elem[5+6]*cy;
  44
  45         anom = elem[6] + elem[6+6]*cy - argp;
  46         motion = elem[6+6] / 36525. / 3600;
  47
  48         incl *= radian;
  49         node *= radian;
  50         argp *= radian;
  51         anom = fmod(anom,360.)*radian;
  52
  53         enom = anom + ecc*sin(anom);
  54         do {
  55                 dele = (anom - enom + ecc * sin(enom)) /
  56                         (1. - ecc*cos(enom));
  57                 enom += dele;
  58         } while(fabs(dele) > converge);
  59         vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  60                 cos(enom/2.));
  61         rad = mrad*(1. - ecc*cos(enom));
  62
  63         lambda = vnom + argp;
  64         pturbl = 0.;
  65         lambda += pturbl*radsec;
  66         pturbb = 0.;
  67         pturbr = 0.;
  68
  69 /*
  70  *      reduce to the ecliptic
  71  */
  72
  73         nd = lambda - node;
  74         lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
  75
  76         sl = sin(incl)*sin(nd) + pturbb*radsec;
  77         beta = atan2(sl, pyth(sl));
  78
  79         lograd = pturbr*2.30258509;
  80         rad *= 1. + lograd;
  81
  82
  83         lambda -= 1185.*radsec;
  84         beta -= 51.*radsec;
  85
  86         motion *= radian*mrad*mrad/(rad*rad);
  87         semi = 83.33;
  88
  89 /*
  90  *      here begins the computation of magnitude
  91  *      first find the geocentric equatorial coordinates of Saturn
  92  */
  93
  94         sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
  95                 sin(beta)*cos(obliq));
  96         sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
  97                 sin(beta)*sin(obliq));
  98         ca = rad*cos(beta)*cos(lambda);
  99         sd += zms;
 100         sa += yms;
 101         ca += xms;
 102         alpha = atan2(sa,ca);
 103         delta = atan2(sd,sqrt(sa*sa+ca*ca));
 104
 105 /*
 106  *      here are the necessary elements of Saturn's rings
 107  *      cf. Exp. Supp. p. 363ff.
 108  */
 109
 110         capj = 6.9056 - 0.4322*capt;
 111         capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
 112         eye = 28.0743 - 0.0128*capt;
 113         comg = 168.1179 + 1.3936*capt;
 114         omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
 115
 116         capj *= radian;
 117         capn *= radian;
 118         eye *= radian;
 119         comg *= radian;
 120         omg *= radian;
 121
 122 /*
 123  *      now find saturnicentric ring-plane coords of the earth
 124  */
 125
 126         sb = sin(capj)*cos(delta)*sin(alpha-capn) -
 127                 cos(capj)*sin(delta);
 128         su = cos(capj)*cos(delta)*sin(alpha-capn) +
 129                 sin(capj)*sin(delta);
 130         cu = cos(delta)*cos(alpha-capn);
 131         u = atan2(su,cu);
 132         b = atan2(sb,sqrt(su*su+cu*cu));
 133
 134 /*
 135  *      and then the saturnicentric ring-plane coords of the sun
 136  */
 137
 138         su = sin(eye)*sin(beta) +
 139                 cos(eye)*cos(beta)*sin(lambda-comg);
 140         cu = cos(beta)*cos(lambda-comg);
 141         up = atan2(su,cu);
 142
 143 /*
 144  *      at last, the magnitude
 145  */
 146
 147
 148         sb = sin(b);
 149         mag = -8.68 +2.52*fabs(up+omg-u)-
 150                 2.60*fabs(sb) + 1.25*(sb*sb);
 151
 152         helio();
 153         geo();
 154 }