URI:
       tmerc.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
       ---
       tmerc.c (1873B)
       ---
            1 #include "astro.h"
            2 
            3 void
            4 merc(void)
            5 {
            6         double pturbl, pturbr;
            7         double lograd;
            8         double dele, enom, vnom, nd, sl;
            9         double q0, v0, t0, j0 , s0;
           10         double lsun, elong, ci, dlong;
           11 
           12         ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;
           13         incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;
           14         node = 47.145944 + 1.185208*capt + .0001739*capt2;
           15         argp = 75.899697 + 1.555490*capt + .0002947*capt2;
           16         mrad = .3870986;
           17         anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;
           18         motion = 4.0923770233;
           19 
           20         q0 = 102.28  + 4.092334429*eday;
           21         v0 = 212.536 + 1.602126105*eday;
           22         t0 = -1.45  + .985604737*eday;
           23         j0 = 225.36 + .083086735*eday;
           24         s0 = 175.68 + .033455441*eday;
           25 
           26         q0 *= radian;
           27         v0 *= radian;
           28         t0 *= radian;
           29         j0 *= radian;
           30         s0 *= radian;
           31 
           32         incl *= radian;
           33         node *= radian;
           34         argp *= radian;
           35         anom = fmod(anom, 360.)*radian;
           36 
           37 
           38         enom = anom + ecc*sin(anom);
           39         do {
           40                 dele = (anom - enom + ecc * sin(enom)) /
           41                         (1. - ecc*cos(enom));
           42                 enom += dele;
           43         } while(fabs(dele) > converge);
           44         vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
           45                         cos(enom/2.));
           46         rad = mrad*(1. - ecc*cos(enom));
           47 
           48         icosadd(mercfp, merccp);
           49         pturbl =  cosadd(2, q0, -v0);
           50         pturbl += cosadd(2, q0, -t0);
           51         pturbl += cosadd(2, q0, -j0);
           52         pturbl += cosadd(2, q0, -s0);
           53 
           54         pturbr =  cosadd(2, q0, -v0);
           55         pturbr += cosadd(2, q0, -t0);
           56         pturbr += cosadd(2, q0, -j0);
           57 
           58 /*
           59  *        reduce to the ecliptic
           60  */
           61 
           62         lambda = vnom + argp + pturbl*radsec;
           63         nd = lambda - node;
           64         lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
           65 
           66         sl = sin(incl)*sin(nd);
           67         beta = atan2(sl, pyth(sl));
           68 
           69         lograd = pturbr*2.30258509;
           70         rad *= 1. + lograd;
           71 
           72         motion *= radian*mrad*mrad/(rad*rad);
           73         semi = 3.34;
           74 
           75         lsun = 99.696678 + 0.9856473354*eday;
           76         lsun *= radian;
           77         elong = lambda - lsun;
           78         ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
           79         dlong = atan2(pyth(ci), ci)/radian;
           80         mag = -.003 + .01815*dlong + .0001023*dlong*dlong;
           81 
           82         helio();
           83         geo();
           84 }