sys/src/cmd/astro/venus.c

   1 #include "astro.h"
   2
   3
   4 void
   5 venus(void)
   6 {
   7         double pturbl, pturbb, pturbr;
   8         double lograd;
   9         double dele, enom, vnom, nd, sl;
  10         double v0, t0, m0, j0, s0;
  11         double lsun, elong, ci, dlong;
  12
  13 /*
  14  *      here are the mean orbital elements
  15  */
  16
  17         ecc = .00682069 - .00004774*capt + 0.091e-6*capt2;
  18         incl = 3.393631 + .0010058*capt - 0.97e-6*capt2;
  19         node = 75.779647 + .89985*capt + .00041*capt2;
  20         argp = 130.163833 + 1.408036*capt - .0009763*capt2;
  21         mrad = .7233316;
  22         anom = 212.603219 + 1.6021301540*eday + .00128605*capt2;
  23         motion = 1.6021687039;
  24
  25 /*
  26  *      mean anomalies of perturbing planets
  27  */
  28
  29         v0 = 212.60 + 1.602130154*eday;
  30         t0 = 358.63  + .985608747*eday;
  31         m0 = 319.74 + 0.524032490*eday;
  32         j0 = 225.43 + .083090842*eday;
  33         s0 = 175.8  + .033459258*eday;
  34
  35         v0 *= radian;
  36         t0 *= radian;
  37         m0 *= radian;
  38         j0 *= radian;
  39         s0 *= radian;
  40
  41         incl *= radian;
  42         node *= radian;
  43         argp *= radian;
  44         anom = fmod(anom, 360.)*radian;
  45
  46 /*
  47  *      computation of long period terms affecting the mean anomaly
  48  */
  49
  50         anom +=
  51                    (2.761-0.022*capt)*radsec*sin(
  52                   13.*t0 - 8.*v0 + 43.83*radian + 4.52*radian*capt)
  53                  + 0.268*radsec*cos(4.*m0 - 7.*t0 + 3.*v0)
  54                  + 0.019*radsec*sin(4.*m0 - 7.*t0 + 3.*v0)
  55                  - 0.208*radsec*sin(s0 + 1.4*radian*capt);
  56
  57 /*
  58  *      computation of elliptic orbit
  59  */
  60
  61         enom = anom + ecc*sin(anom);
  62         do {
  63                 dele = (anom - enom + ecc * sin(enom)) /
  64                         (1 - ecc*cos(enom));
  65                 enom += dele;
  66         } while(fabs(dele) > converge);
  67         vnom = 2*atan2(sqrt((1+ecc)/(1-ecc))*sin(enom/2),
  68                 cos(enom/2));
  69         rad = mrad*(1 - ecc*cos(enom));
  70
  71         lambda = vnom + argp;
  72
  73 /*
  74  *      perturbations in longitude
  75  */
  76
  77         icosadd(venfp, vencp);
  78         pturbl = cosadd(4, v0, t0, m0, j0);
  79         pturbl *= radsec;
  80
  81 /*
  82  *      perturbations in latidude
  83  */
  84
  85         pturbb = cosadd(3, v0, t0, j0);
  86         pturbb *= radsec;
  87
  88 /*
  89  *      perturbations in log radius vector
  90  */
  91
  92         pturbr = cosadd(4, v0, t0, m0, j0);
  93
  94 /*
  95  *      reduction to the ecliptic
  96  */
  97
  98         lambda += pturbl;
  99         nd = lambda - node;
 100         lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
 101
 102         sl = sin(incl)*sin(nd);
 103         beta = atan2(sl, pyth(sl)) + pturbb;
 104
 105         lograd = pturbr*2.30258509;
 106         rad *= 1 + lograd;
 107
 108
 109         motion *= radian*mrad*mrad/(rad*rad);
 110
 111 /*
 112  *      computation of magnitude
 113  */
 114
 115         lsun = 99.696678 + 0.9856473354*eday;
 116         lsun *= radian;
 117         elong = lambda - lsun;
 118         ci = (rad - cos(elong))/sqrt(1 + rad*rad - 2*rad*cos(elong));
 119         dlong = atan2(pyth(ci), ci)/radian;
 120         mag = -4 + .01322*dlong + .0000004247*dlong*dlong*dlong;
 121
 122         semi = 8.41;
 123
 124         helio();
 125         geo();
 126 }