sys/src/cmd/astro/mars.c

   1 #include "astro.h"
   2
   3 void
   4 mars(void)
   5 {
   6         double pturbl, pturbb, pturbr;
   7         double lograd;
   8         double dele, enom, vnom, nd, sl;
   9         double lsun, elong, ci, dlong;
  10
  11
  12         ecc = .09331290 + .000092064*capt;
  13         incl = 1.850333 - 6.75e-4*capt;
  14         node = 48.786442 + .770992*capt;
  15         argp = 334.218203 + 1.840758*capt + 1.30e-4*capt2;
  16         mrad = 1.5236915;
  17         anom = 319.529425 + .5240207666*eday + 1.808e-4*capt2;
  18         motion = 0.5240711638;
  19
  20
  21         incl = incl*radian;
  22         node = node*radian;
  23         argp = argp*radian;
  24         anom = fmod(anom,360.)*radian;
  25
  26         enom = anom + ecc*sin(anom);
  27         do {
  28                 dele = (anom - enom + ecc * sin(enom)) /
  29                         (1. - ecc*cos(enom));
  30                 enom += dele;
  31         } while(fabs(dele) > converge);
  32         vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),cos(enom/2.));
  33         rad = mrad*(1. - ecc*cos(enom));
  34
  35         lambda = vnom + argp;
  36         pturbl = 0.;
  37         lambda = lambda + pturbl*radsec;
  38         pturbb = 0.;
  39         pturbr = 0.;
  40
  41 /*
  42  *      reduce to the ecliptic
  43  */
  44
  45         nd = lambda - node;
  46         lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
  47
  48         sl = sin(incl)*sin(nd) + pturbb*radsec;
  49         beta = atan2(sl, pyth(sl));
  50
  51         lograd = pturbr*2.30258509;
  52         rad *= 1. + lograd;
  53
  54
  55         motion *= radian*mrad*mrad/(rad*rad);
  56         semi = 4.68;
  57
  58         lsun = 99.696678 + 0.9856473354*eday;
  59         lsun *= radian;
  60         elong = lambda - lsun;
  61         ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
  62         dlong = atan2(pyth(ci), ci)/radian;
  63         mag = -1.30 + .01486*dlong;
  64
  65         helio();
  66         geo();
  67 }