tvenus.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tvenus.c (2417B)
       ---
            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 }