home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume9 / ephem2 / part04 < prev    next >
Encoding:
Text File  |  1989-11-27  |  44.8 KB  |  1,629 lines
  1. Newsgroups: comp.sources.misc
  2. From: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  3. Subject: v09i041: ephem, v4.8, 4 of 5
  4. Reply-To: ecd@umn-cs.cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
  5.  
  6. Posting-number: Volume 9, Issue 41
  7. Submitted-by: ecd@umn-cs.cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
  8. Archive-name: ephem2/part04
  9.  
  10. # This is a "shell archive" file; run it with sh.
  11. # This is file 4.
  12. echo x plans.c
  13. cat > plans.c << 'xXx'
  14. #include <stdio.h>
  15. #include <math.h>
  16. #include "astro.h"
  17.  
  18. #define    TWOPI        (2*PI)
  19. #define    mod2PI(x)    ((x) - (int)((x)/TWOPI)*TWOPI)
  20.  
  21. /* given a modified Julian date, mjd, and a planet, p, find:
  22.  *   lpd0: heliocentric longitude, 
  23.  *   psi0: heliocentric latitude,
  24.  *   rp0:  distance from the sun to the planet, 
  25.  *   rho0: distance from the Earth to the planet,
  26.  *         none corrected for light time, ie, they are the true values for the
  27.  *         given instant.
  28.  *   lam:  geocentric ecliptic longitude, 
  29.  *   bet:  geocentric ecliptic latitude,
  30.  *         each corrected for light time, ie, they are the apparent values as
  31.  *       seen from the center of the Earth for the given instant.
  32.  *   dia:  angular diameter in arcsec at 1 AU, 
  33.  *   mag:  visual magnitude when 1 AU from sun and earth at 0 phase angle.
  34.  *
  35.  * all angles are in radians, all distances in AU.
  36.  * the mean orbital elements are found by calling pelement(), then mutual
  37.  *   perturbation corrections are applied as necessary.
  38.  *
  39.  * corrections for nutation and abberation must be made by the caller. The RA 
  40.  *   and DEC calculated from the fully-corrected ecliptic coordinates are then
  41.  *   the apparent geocentric coordinates. Further corrections can be made, if
  42.  *   required, for atmospheric refraction and geocentric parallax although the
  43.  *   intrinsic error herein of about 10 arcseconds is usually the dominant
  44.  *   error at this stage.
  45.  * TODO: combine the several intermediate expressions when get a good compiler.
  46.  */
  47. plans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
  48. double mjd;
  49. int p;
  50. double *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
  51. {
  52.     static double plan[8][9];
  53.     static lastmjd = -10000;
  54.     double dl;    /* perturbation correction for longitude */
  55.     double dr;    /*  "   orbital radius */
  56.     double dml;    /*  "   mean longitude */
  57.     double ds;    /*  "   eccentricity */
  58.     double dm;    /*  "   mean anomaly */
  59.     double da;    /*  "   semi-major axis */
  60.     double dhl;    /*  "   heliocentric longitude */
  61.     double lsn, rsn;    /* true geocentric longitude of sun and sun-earth rad */
  62.     double mas;    /* mean anomaly of the sun */
  63.     double re;    /* radius of earth's orbit */
  64.     double lg;    /* longitude of earth */
  65.     double map[8];    /* array of mean anomalies for each planet */
  66.     double lpd, psi, rp, rho;
  67.     double ll, sll, cll;
  68.     double t;
  69.     double dt;
  70.     int pass;
  71.     int j;
  72.     double s, ma;
  73.     double nu, ea;
  74.     double lp, om;
  75.     double lo, slo, clo;
  76.     double inc, y;
  77.     double spsi, cpsi;
  78.     double rpd;
  79.  
  80.     /* only need to fill in plan[] once for a given mjd */
  81.     if (mjd != lastmjd) {
  82.         pelement (mjd, plan);
  83.         lastmjd = mjd;
  84.     }
  85.  
  86.     dt = 0;
  87.     t = mjd/36525.;
  88.     sunpos (mjd, &lsn, &rsn);
  89.     masun (mjd, &mas);
  90.         re = rsn;
  91.     lg = lsn+PI;
  92.  
  93.     /* first find the true position of the planet at mjd.
  94.      * then repeat a second time for a slightly different time based
  95.      * on the position found in the first pass to account for light-travel
  96.      * time.
  97.      */
  98.     for (pass = 0; pass < 2; pass++) {
  99.  
  100.         for (j = 0; j < 8; j++)
  101.         map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]);
  102.  
  103.         /* set initial corrections to 0.
  104.          * then modify as necessary for the planet of interest.
  105.          */
  106.         dl = 0;
  107.         dr = 0;
  108.         dml = 0;
  109.         ds = 0;
  110.         dm = 0;
  111.         da = 0;
  112.         dhl = 0;
  113.  
  114.         switch (p) {
  115.  
  116.         case MERCURY:
  117.         p_mercury (map, &dl, &dr);
  118.         break;
  119.  
  120.         case VENUS:
  121.         p_venus (t, mas, map, &dl, &dr, &dml, &dm);
  122.         break;
  123.  
  124.         case MARS:
  125.         p_mars (mas, map, &dl, &dr, &dml, &dm);
  126.         break;
  127.  
  128.         case JUPITER:
  129.         p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da);
  130.         break;
  131.  
  132.         case SATURN:
  133.         p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl);
  134.         break;
  135.  
  136.         case URANUS:
  137.         p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
  138.         break;
  139.  
  140.         case NEPTUNE:
  141.         p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
  142.         break;
  143.  
  144.         case PLUTO:
  145.         /* no perturbation theory for pluto */
  146.         break;
  147.         }
  148.  
  149.         s = plan[p][3]+ds;
  150.         ma = map[p]+dm;
  151.         anomaly (ma, s, &nu, &ea);
  152.         rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu));
  153.         lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm);
  154.         lp = degrad(lp);
  155.         om = degrad(plan[p][5]);
  156.         lo = lp-om;
  157.         slo = sin(lo);
  158.         clo = cos(lo);
  159.         inc = degrad(plan[p][4]);
  160.         rp = rp+dr;
  161.         spsi = slo*sin(inc);
  162.         y = slo*cos(inc);
  163.         psi = asin(spsi)+dhl;
  164.         spsi = sin(psi);
  165.         lpd = atan(y/clo)+om+degrad(dl);
  166.         if (clo<0) lpd += PI;
  167.         if (lpd>TWOPI) lpd -= TWOPI;
  168.         cpsi = cos(psi);
  169.         rpd = rp*cpsi;
  170.         ll = lpd-lg;
  171.         rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll));
  172.  
  173.         /* when we view a planet we see it in the position it occupied
  174.          * dt hours ago, where rho is the distance between it and earth,
  175.          * in AU. use this as the new time for the next pass.
  176.          */
  177.         dt = rho*5.775518e-3;
  178.  
  179.         if (pass == 0) {
  180.         /* save heliocentric coordinates after first pass since, being
  181.          * true, they are NOT to be corrected for light-travel time.
  182.          */
  183.         *lpd0 = lpd;
  184.         *psi0 = psi;
  185.         *rp0 = rp;
  186.         *rho0 = rho;
  187.         }
  188.     }
  189.  
  190.         sll = sin(ll);
  191.     cll = cos(ll);
  192.         if (p < MARS) 
  193.         *lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI;
  194.     else
  195.         *lam = atan(re*sll/(rpd-re*cll))+lpd;
  196.     range (lam, TWOPI);
  197.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll));
  198.     *dia = plan[p][7];
  199.     *mag = plan[p][8];
  200. }
  201.  
  202. /* set auxilliary variables used for jupiter, saturn, uranus, and neptune */
  203. static
  204. aux_jsun (t, j1, j2, j3, j4, j5, j6)
  205. double t;
  206. double *j1, *j2, *j3, *j4, *j5, *j6;
  207. {
  208.         *j1 = t/5+0.1;
  209.         *j2 = mod2PI(4.14473+5.29691e1*t);
  210.         *j3 = mod2PI(4.641118+2.132991e1*t);
  211.         *j4 = mod2PI(4.250177+7.478172*t);
  212.         *j5 = 5 * *j3 - 2 * *j2;
  213.     *j6 = 2 * *j2 - 6 * *j3 + 3 * *j4;
  214. }
  215.  
  216. /* find the mean anomaly of the sun at mjd.
  217.  * this is the same as that used in sun() but when it was converted to C it
  218.  * was not known it would be required outside that routine.
  219.  * TODO: add an argument to sun() to return mas and eliminate this routine.
  220.  */
  221. static
  222. masun (mjd, mas)
  223. double mjd;
  224. double *mas;
  225. {
  226.     double t, t2;
  227.     double a, b;
  228.  
  229.     t = mjd/36525;
  230.     t2 = t*t;
  231.     a = 9.999736042e1*t;
  232.     b = 360.*(a-(int)a);
  233.     *mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b);
  234. }
  235.  
  236. /* perturbations for mercury */
  237. static
  238. p_mercury (map, dl, dr)
  239. double map[];
  240. double *dl, *dr;
  241. {
  242.     *dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+
  243.          1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+
  244.          9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+
  245.          7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1);
  246.  
  247.     *dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+
  248.          6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+
  249.          5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+
  250.          3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699);
  251. }
  252.  
  253. /* ....venus */
  254. static
  255. p_venus (t, mas, map, dl, dr, dml, dm)
  256. double t, mas, map[];
  257. double *dl, *dr, *dml, *dm;
  258. {
  259.     *dml = degrad (7.7e-4*sin(4.1406+t*2.6227));
  260.     *dm = *dml;
  261.  
  262.     *dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+
  263.          1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+
  264.          1.36e-3*cos(mas-map[2-1]-2.0788)+
  265.          9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+
  266.          8.2e-4*cos(map[3]-map[2-1]-3.6318);
  267.  
  268.     *dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+
  269.          1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+
  270.          6.887e-6*cos(map[3]-map[2-1]-2.06106)+
  271.          5.172e-6*cos(mas-map[2-1]-5.08065e-1)+
  272.          3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+
  273.          3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+
  274.          3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1);
  275. }
  276.  
  277. /* ....mars */
  278. static
  279. p_mars (mas, map, dl, dr, dml, dm)
  280. double mas, map[];
  281. double *dl, *dr, *dml, *dm;
  282. {
  283.     double a;
  284.  
  285.     a = 3*map[3]-8*map[2]+4*mas;
  286.     *dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
  287.     *dm = *dml;
  288.  
  289.     *dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
  290.          6.07e-3*cos(2*map[3]-map[2]-3.2873)+
  291.          4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
  292.          3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
  293.          2.38e-3*cos(mas-map[2]+6.1256e-1)+
  294.          2.04e-3*cos(2*mas-3*map[2]+2.7688)+
  295.          1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
  296.          1.36e-3*cos(2*mas-4*map[2]+2.6894)+
  297.          1.04e-3*cos(map[3]+3.0749e-1);
  298.  
  299.     *dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
  300.          5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
  301.          3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
  302.          1.5996e-5*cos(mas-map[2]-9.69618e-1)+
  303.          1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
  304.          8.966e-6*cos(map[3]-2*map[2]+7.61225e-1);
  305.      *dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
  306.          7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
  307.          6.62e-6*cos(mas-2*map[2]+1.97575)+
  308.          4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
  309.          4.693e-6*cos(3*mas-5*map[2]+3.32665)+
  310.          4.571e-6*cos(2*mas-4*map[2]+4.27086)+
  311.          4.409e-6*cos(3*map[3]-map[2]-2.02158);
  312. }
  313.  
  314. /* ....jupiter */
  315. static
  316. p_jupiter (t, s, dml, ds, dm, da)
  317. double t, s;
  318. double *dml, *ds, *dm, *da;
  319. {
  320.     double dp;
  321.     double j1, j2, j3, j4, j5, j6, j7;
  322.     double sj3, cj3, s2j3, c2j3;
  323.         double sj5, cj5, s2j5;
  324.     double sj6;
  325.         double sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7;
  326.  
  327.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  328.         j7 = j3-j2;
  329.     sj3 = sin(j3);
  330.     cj3 = cos(j3);
  331.         s2j3 = sin(2*j3);
  332.     c2j3 = cos(2*j3);
  333.         sj5 = sin(j5);
  334.     cj5 = cos(j5);
  335.         s2j5 = sin(2*j5);
  336.     sj6 = sin(j6);
  337.         sj7 = sin(j7);
  338.     cj7 = cos(j7);
  339.         s2j7 = sin(2*j7);
  340.     c2j7 = cos(2*j7);
  341.         s3j7 = sin(3*j7);
  342.     c3j7 = cos(3*j7);
  343.         s4j7 = sin(4*j7);
  344.     c4j7 = cos(4*j7);
  345.         c5j7 = cos(5*j7);
  346.  
  347.     *dml = (3.31364e-1-(1.0281e-2+4.692e-3*j1)*j1)*sj5+
  348.           (3.228e-3-(6.4436e-2-2.075e-3*j1)*j1)*cj5-
  349.           (3.083e-3+(2.75e-4-4.89e-4*j1)*j1)*s2j5+
  350.           2.472e-3*sj6+1.3619e-2*sj7+1.8472e-2*s2j7+6.717e-3*s3j7+
  351.           2.775e-3*s4j7+6.417e-3*s2j7*sj3+
  352.           (7.275e-3-1.253e-3*j1)*sj7*sj3+
  353.           2.439e-3*s3j7*sj3-(3.5681e-2+1.208e-3*j1)*sj7*cj3;
  354.         *dml += -3.767e-3*c2j7*sj3-(3.3839e-2+1.125e-3*j1)*cj7*sj3-
  355.           4.261e-3*s2j7*cj3+
  356.           (1.161e-3*j1-6.333e-3)*cj7*cj3+
  357.           2.178e-3*cj3-6.675e-3*c2j7*cj3-2.664e-3*c3j7*cj3-
  358.           2.572e-3*sj7*s2j3-3.567e-3*s2j7*s2j3+2.094e-3*cj7*c2j3+
  359.           3.342e-3*c2j7*c2j3;
  360.     *dml = degrad(*dml);
  361.  
  362.     *ds = (3606+(130-43*j1)*j1)*sj5+(1289-580*j1)*cj5-6764*sj7*sj3-
  363.          1110*s2j7*sj3-224*s3j7*sj3-204*sj3+(1284+116*j1)*cj7*sj3+
  364.          188*c2j7*sj3+(1460+130*j1)*sj7*cj3+224*s2j7*cj3-817*cj3+
  365.          6074*cj3*cj7+992*c2j7*cj3+
  366.          508*c3j7*cj3+230*c4j7*cj3+108*c5j7*cj3;
  367.     *ds += -(956+73*j1)*sj7*s2j3+448*s2j7*s2j3+137*s3j7*s2j3+
  368.          (108*j1-997)*cj7*s2j3+480*c2j7*s2j3+148*c3j7*s2j3+
  369.          (99*j1-956)*sj7*c2j3+490*s2j7*c2j3+
  370.          158*s3j7*c2j3+179*c2j3+(1024+75*j1)*cj7*c2j3-
  371.          437*c2j7*c2j3-132*c3j7*c2j3;
  372.     *ds *= 1e-7;
  373.  
  374.     dp = (7.192e-3-3.147e-3*j1)*sj5-4.344e-3*sj3+
  375.          (j1*(1.97e-4*j1-6.75e-4)-2.0428e-2)*cj5+
  376.          3.4036e-2*cj7*sj3+(7.269e-3+6.72e-4*j1)*sj7*sj3+
  377.          5.614e-3*c2j7*sj3+2.964e-3*c3j7*sj3+3.7761e-2*sj7*cj3+
  378.          6.158e-3*s2j7*cj3-
  379.          6.603e-3*cj7*cj3-5.356e-3*sj7*s2j3+2.722e-3*s2j7*s2j3+
  380.          4.483e-3*cj7*s2j3-2.642e-3*c2j7*s2j3+4.403e-3*sj7*c2j3-
  381.          2.536e-3*s2j7*c2j3+5.547e-3*cj7*c2j3-2.689e-3*c2j7*c2j3;
  382.  
  383.     *dm = *dml-(degrad(dp)/s);
  384.  
  385.     *da = 205*cj7-263*cj5+693*c2j7+312*c3j7+147*c4j7+299*sj7*sj3+
  386.          181*c2j7*sj3+204*s2j7*cj3+111*s3j7*cj3-337*cj7*cj3-
  387.          111*c2j7*cj3;
  388.     *da *= 1e-6;
  389. }
  390.  
  391. /* ....saturn */
  392. static
  393. p_saturn (t, s, dml, ds, dm, da, dhl)
  394. double t, s;
  395. double *dml, *ds, *dm, *da, *dhl;
  396. {
  397.     double dp;
  398.     double j1, j2, j3, j4, j5, j6, j7, j8;
  399.     double sj3, cj3, s2j3, c2j3, s3j3, c3j3, s4j3, c4j3;
  400.         double sj5, cj5, s2j5, c2j5;
  401.     double sj6;
  402.         double sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7, s5j7;
  403.     double s2j8, c2j8, s3j8, c3j8;
  404.  
  405.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  406.         j7 = j3-j2;
  407.     sj3 = sin(j3);
  408.     cj3 = cos(j3);
  409.         s2j3 = sin(2*j3);
  410.     c2j3 = cos(2*j3);
  411.         sj5 = sin(j5);
  412.     cj5 = cos(j5);
  413.         s2j5 = sin(2*j5);
  414.     sj6 = sin(j6);
  415.         sj7 = sin(j7);
  416.     cj7 = cos(j7);
  417.         s2j7 = sin(2*j7);
  418.     c2j7 = cos(2*j7);
  419.         s3j7 = sin(3*j7);
  420.     c3j7 = cos(3*j7);
  421.         s4j7 = sin(4*j7);
  422.     c4j7 = cos(4*j7);
  423.         c5j7 = cos(5*j7);
  424.  
  425.     s3j3 = sin(3*j3);
  426.     c3j3 = cos(3*j3);
  427.     s4j3 = sin(4*j3);
  428.     c4j3 = cos(4*j3);
  429.     c2j5 = cos(2*j5);
  430.     s5j7 = sin(5*j7);
  431.     j8 = j4-j3;
  432.     s2j8 = sin(2*j8);
  433.     c2j8 = cos(2*j8);
  434.     s3j8 = sin(3*j8);
  435.     c3j8 = cos(3*j8);
  436.  
  437.     *dml = 7.581e-3*s2j5-7.986e-3*sj6-1.48811e-1*sj7-4.0786e-2*s2j7-
  438.           (8.14181e-1-(1.815e-2-1.6714e-2*j1)*j1)*sj5-
  439.           (1.0497e-2-(1.60906e-1-4.1e-3*j1)*j1)*cj5-1.5208e-2*s3j7-
  440.           6.339e-3*s4j7-6.244e-3*sj3-1.65e-2*s2j7*sj3+
  441.           (8.931e-3+2.728e-3*j1)*sj7*sj3-5.775e-3*s3j7*sj3+
  442.           (8.1344e-2+3.206e-3*j1)*cj7*sj3+1.5019e-2*c2j7*sj3;
  443.     *dml += (8.5581e-2+2.494e-3*j1)*sj7*cj3+1.4394e-2*c2j7*cj3+
  444.           (2.5328e-2-3.117e-3*j1)*cj7*cj3+
  445.           6.319e-3*c3j7*cj3+6.369e-3*sj7*s2j3+9.156e-3*s2j7*s2j3+
  446.           7.525e-3*s3j8*s2j3-5.236e-3*cj7*c2j3-7.736e-3*c2j7*c2j3-
  447.           7.528e-3*c3j8*c2j3;
  448.     *dml = degrad(*dml);
  449.  
  450.     *ds = (-7927+(2548+91*j1)*j1)*sj5+(13381+(1226-253*j1)*j1)*cj5+
  451.          (248-121*j1)*s2j5-(305+91*j1)*c2j5+412*s2j7+12415*sj3+
  452.          (390-617*j1)*sj7*sj3+(165-204*j1)*s2j7*sj3+26599*cj7*sj3-
  453.          4687*c2j7*sj3-1870*c3j7*sj3-821*c4j7*sj3-
  454.          377*c5j7*sj3+497*c2j8*sj3+(163-611*j1)*cj3;
  455.     *ds += -12696*sj7*cj3-4200*s2j7*cj3-1503*s3j7*cj3-619*s4j7*cj3-
  456.          268*s5j7*cj3-(282+1306*j1)*cj7*cj3+(-86+230*j1)*c2j7*cj3+
  457.          461*s2j8*cj3-350*s2j3+(2211-286*j1)*sj7*s2j3-
  458.          2208*s2j7*s2j3-568*s3j7*s2j3-346*s4j7*s2j3-
  459.          (2780+222*j1)*cj7*s2j3+(2022+263*j1)*c2j7*s2j3+248*c3j7*s2j3+
  460.          242*s3j8*s2j3+467*c3j8*s2j3-490*c2j3-(2842+279*j1)*sj7*c2j3;
  461.     *ds += (128+226*j1)*s2j7*c2j3+224*s3j7*c2j3+
  462.          (-1594+282*j1)*cj7*c2j3+(2162-207*j1)*c2j7*c2j3+
  463.          561*c3j7*c2j3+343*c4j7*c2j3+469*s3j8*c2j3-242*c3j8*c2j3-
  464.          205*sj7*s3j3+262*s3j7*s3j3+208*cj7*c3j3-271*c3j7*c3j3-
  465.          382*c3j7*s4j3-376*s3j7*c4j3;
  466.     *ds *= 1e-7;
  467.  
  468.     dp = (7.7108e-2+(7.186e-3-1.533e-3*j1)*j1)*sj5-7.075e-3*sj7+
  469.          (4.5803e-2-(1.4766e-2+5.36e-4*j1)*j1)*cj5-7.2586e-2*cj3-
  470.          7.5825e-2*sj7*sj3-2.4839e-2*s2j7*sj3-8.631e-3*s3j7*sj3-
  471.          1.50383e-1*cj7*cj3+2.6897e-2*c2j7*cj3+1.0053e-2*c3j7*cj3-
  472.          (1.3597e-2+1.719e-3*j1)*sj7*s2j3+1.1981e-2*s2j7*c2j3;
  473.     dp += -(7.742e-3-1.517e-3*j1)*cj7*s2j3+
  474.          (1.3586e-2-1.375e-3*j1)*c2j7*c2j3-
  475.          (1.3667e-2-1.239e-3*j1)*sj7*c2j3+
  476.          (1.4861e-2+1.136e-3*j1)*cj7*c2j3-
  477.          (1.3064e-2+1.628e-3*j1)*c2j7*c2j3;
  478.  
  479.     *dm = *dml-(degrad(dp)/s);
  480.  
  481.     *da = 572*sj5-1590*s2j7*cj3+2933*cj5-647*s3j7*cj3+33629*cj7-
  482.          344*s4j7*cj3-3081*c2j7+2885*cj7*cj3-1423*c3j7+
  483.          (2172+102*j1)*c2j7*cj3-671*c4j7+296*c3j7*cj3-320*c5j7-
  484.          267*s2j7*s2j3+1098*sj3-778*cj7*s2j3-2812*sj7*sj3;
  485.     *da += 495*c2j7*s2j3+688*s2j7*sj3+250*c3j7*s2j3-393*s3j7*sj3-
  486.          856*sj7*c2j3-228*s4j7*sj3+441*s2j7*c2j3+2138*cj7*sj3+
  487.          296*c2j7*c2j3-999*c2j7*sj3+211*c3j7*c2j3-642*c3j7*sj3-
  488.          427*sj7*s3j3-325*c4j7*sj3+398*s3j7*s3j3-890*cj3+
  489.          344*cj7*c3j3+2206*sj7*cj3-427*c3j7*c3j3;
  490.     *da *= 1e-6;
  491.  
  492.     *dhl = 7.47e-4*cj7*sj3+1.069e-3*cj7*cj3+2.108e-3*s2j7*s2j3+
  493.           1.261e-3*c2j7*s2j3+1.236e-3*s2j7*c2j3-2.075e-3*c2j7*c2j3;
  494.     *dhl = degrad(*dhl);
  495. }
  496.  
  497. /* ....uranus */
  498. static
  499. p_uranus (t, s, dl, dr, dml, ds, dm, da, dhl)
  500. double t, s;
  501. double *dl, *dr, *dml, *ds, *dm, *da, *dhl;
  502. {
  503.     double dp;
  504.     double j1, j2, j3, j4, j5, j6;
  505.     double j8, j9, j10, j11, j12;
  506.     double sj4, cj4, s2j4, c2j4;
  507.     double sj9, cj9, s2j9, c2j9;
  508.     double sj11, cj11;
  509.  
  510.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  511.  
  512.         j8 = mod2PI(1.46205+3.81337*t);
  513.         j9 = 2*j8-j4;
  514.     sj9 = sin(j9);
  515.     cj9 = cos(j9);
  516.         s2j9 = sin(2*j9);
  517.     c2j9 = cos(2*j9);
  518.  
  519.     j10 = j4-j2;
  520.     j11 = j4-j3;
  521.     j12 = j8-j4;
  522.  
  523.     *dml = (8.64319e-1-1.583e-3*j1)*sj9+(8.2222e-2-6.833e-3*j1)*cj9+
  524.           3.6017e-2*s2j9-3.019e-3*c2j9+8.122e-3*sin(j6);
  525.     *dml = degrad(*dml);
  526.  
  527.     dp = 1.20303e-1*sj9+6.197e-3*s2j9+(1.9472e-2-9.47e-4*j1)*cj9;
  528.     *dm = *dml-(degrad(dp)/s);
  529.  
  530.     *ds = (163*j1-3349)*sj9+20981*cj9+1311*c2j9;
  531.     *ds *= 1e-7;
  532.  
  533.     *da = -3.825e-3*cj9;
  534.  
  535.     *dl = (1.0122e-2-9.88e-4*j1)*sin(j4+j11)+
  536.          (-3.8581e-2+(2.031e-3-1.91e-3*j1)*j1)*cos(j4+j11)+
  537.          (3.4964e-2-(1.038e-3-8.68e-4*j1)*j1)*cos(2*j4+j11)+
  538.          5.594e-3*sin(j4+3*j12)-1.4808e-2*sin(j10)-
  539.          5.794e-3*sin(j11)+2.347e-3*cos(j11)+9.872e-3*sin(j12)+
  540.          8.803e-3*sin(2*j12)-4.308e-3*sin(3*j12);
  541.  
  542.     sj11 = sin(j11);
  543.     cj11 = cos(j11);
  544.     sj4 = sin(j4);
  545.     cj4 = cos(j4);
  546.     s2j4 = sin(2*j4);
  547.     c2j4 = cos(2*j4);
  548.     *dhl = (4.58e-4*sj11-6.42e-4*cj11-5.17e-4*cos(4*j12))*sj4-
  549.           (3.47e-4*sj11+8.53e-4*cj11+5.17e-4*sin(4*j11))*cj4+
  550.           4.03e-4*(cos(2*j12)*s2j4+sin(2*j12)*c2j4);
  551.     *dhl = degrad(*dhl);
  552.  
  553.     *dr = -25948+4985*cos(j10)-1230*cj4+3354*cos(j11)+904*cos(2*j12)+
  554.          894*(cos(j12)-cos(3*j12))+(5795*cj4-1165*sj4+1388*c2j4)*sj11+
  555.          (1351*cj4+5702*sj4+1388*s2j4)*cos(j11);
  556.     *dr *= 1e-6;
  557. }
  558.  
  559. /* ....neptune */
  560. static
  561. p_neptune (t, s, dl, dr, dml, ds, dm, da, dhl)
  562. double t, s;
  563. double *dl, *dr, *dml, *ds, *dm, *da, *dhl;
  564. {
  565.     double dp;
  566.     double j1, j2, j3, j4, j5, j6;
  567.     double j8, j9, j10, j11, j12;
  568.     double sj8, cj8;
  569.     double sj9, cj9, s2j9, c2j9;
  570.     double s2j12, c2j12;
  571.  
  572.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  573.  
  574.         j8 = mod2PI(1.46205+3.81337*t);
  575.         j9 = 2*j8-j4;
  576.     sj9 = sin(j9);
  577.     cj9 = cos(j9);
  578.         s2j9 = sin(2*j9);
  579.     c2j9 = cos(2*j9);
  580.  
  581.     j10 = j8-j2;
  582.     j11 = j8-j3;
  583.     j12 = j8-j4;
  584.  
  585.     *dml = (1.089e-3*j1-5.89833e-1)*sj9+(4.658e-3*j1-5.6094e-2)*cj9-
  586.           2.4286e-2*s2j9;
  587.     *dml = degrad(*dml);
  588.  
  589.     dp = 2.4039e-2*sj9-2.5303e-2*cj9+6.206e-3*s2j9-5.992e-3*c2j9;
  590.  
  591.     *dm = *dml-(degrad(dp)/s);
  592.  
  593.     *ds = 4389*sj9+1129*s2j9+4262*cj9+1089*c2j9;
  594.     *ds *= 1e-7;
  595.  
  596.     *da = 8189*cj9-817*sj9+781*c2j9;
  597.     *da *= 1e-6;
  598.  
  599.     s2j12 = sin(2*j12);
  600.     c2j12 = cos(2*j12);
  601.     sj8 = sin(j8);
  602.     cj8 = cos(j8);
  603.     *dl = -9.556e-3*sin(j10)-5.178e-3*sin(j11)+2.572e-3*s2j12-
  604.          2.972e-3*c2j12*sj8-2.833e-3*s2j12*cj8;
  605.  
  606.     *dhl = 3.36e-4*c2j12*sj8+3.64e-4*s2j12*cj8;
  607.     *dhl = degrad(*dhl);
  608.  
  609.     *dr = -40596+4992*cos(j10)+2744*cos(j11)+2044*cos(j12)+1051*c2j12;
  610.     *dr *= 1e-6;
  611. }
  612. xXx
  613. echo x plot.c
  614. cat > plot.c << 'xXx'
  615. /* code to support the plotting capabilities.
  616.  * idea is to let the operator name a plot file and mark some fields for
  617.  * logging. then after each screen update, the logged fields are written to
  618.  * the plot file. later, the file may be plotted (very simplistically by 
  619.  * ephem, for now anyway, or by some other program entirely.).
  620.  * 
  621.  * format of the plot file is one line per coordinate: label,x,y
  622.  * if z was specified, it is a fourth field.
  623.  * x,y,z are plotted using %g format.
  624.  */
  625.  
  626. #include <stdio.h>
  627. #include <math.h>
  628. #include "screen.h"
  629.  
  630. extern errno;
  631. extern char *sys_errlist[];
  632. #define    errsys    (sys_errlist[errno])
  633.  
  634. #define    TRACE(x)    {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
  635.  
  636. #define    MAXPLTLINES    10    /* max number of labeled lines we can track.
  637.                  * note we can't store more than NFLOGS fields
  638.                  * anyway (see flog.c).
  639.                  */
  640. #define    FNLEN        (14+1)    /* longest filename; plus 1 for \0 */
  641.  
  642. static char plt_filename[FNLEN] = "ephem.plt";    /* default plot file name */
  643. static FILE *plt_fp;        /* the plot file; == 0 means don't plot */
  644.  
  645. /* store the label and rcfpack()s for each line to track. */
  646. typedef struct {
  647.     char pl_label;
  648.     int pl_rcpx, pl_rcpy, pl_rcpz;
  649. } PltLine;
  650. static PltLine pltlines[MAXPLTLINES];
  651. static int npltlines;        /* number of pltlines[] in actual use */
  652.  
  653. static int plt_in_polar;    /*if true plot in polar coords, else cartesian*/
  654.  
  655. /* picked the Plot label:
  656.  * if on, just turn it off.
  657.  * if off, turn on, define fields or select name of file to plot and do it.
  658.  * TODO: more flexability, more relevance.
  659.  */
  660. plot_setup()
  661. {
  662.     if (plt_fp) {
  663.         plt_turn_off();
  664.         plot_prstate (0);
  665.     } else {
  666.         static char *chcs[4] = {
  667.         "Select fields", "Display a plot file", (char *)0,
  668.         "Begin plotting"
  669.         };
  670.         int ff = 0;
  671.     ask:
  672.         chcs[2] = plt_in_polar ? "Polar coords" : "Cartesian coords";
  673.         switch (popup(chcs, ff, npltlines > 0 ? 4 : 3)) {
  674.         case 0: plt_select_fields(); break;
  675.         case 1: plt_file(); break;
  676.         case 2: plt_in_polar ^= 1; ff = 2; goto ask;
  677.         case 3: plt_turn_on(); break;
  678.         }
  679.     }
  680. }
  681.  
  682. /* write the active plotfields to the current plot file, if one is open. */
  683. plot()
  684. {
  685.     if (plt_fp) {
  686.         PltLine *plp;
  687.         double x, y, z;
  688.         for (plp = &pltlines[npltlines]; --plp >= pltlines; ) {
  689.         if (flog_get (plp->pl_rcpx, &x) == 0 
  690.             && flog_get (plp->pl_rcpy, &y) == 0) {
  691.             fprintf (plt_fp, "%c,%.12g,%.12g", plp->pl_label, x, y);
  692.             if (flog_get (plp->pl_rcpz, &z) == 0)
  693.             fprintf (plt_fp, ",%.12g", z);
  694.             fprintf (plt_fp, "\n");
  695.         }
  696.         }
  697.     }
  698. }
  699.  
  700. plot_prstate (force)
  701. int force;
  702. {
  703.     static last;
  704.     int this = plt_fp != 0;
  705.  
  706.     if (force || this != last) {
  707.         f_string (R_PLOT, C_PLOTV, this ? " on" : "off");
  708.         last = this;
  709.     }
  710. }
  711.  
  712. plot_ison()
  713. {
  714.     return (plt_fp != 0);
  715. }
  716.  
  717. static
  718. plt_reset()
  719. {
  720.     PltLine *plp;
  721.  
  722.     for (plp = &pltlines[npltlines]; --plp >= pltlines; ) {
  723.         (void) flog_delete (plp->pl_rcpx);
  724.         (void) flog_delete (plp->pl_rcpy);
  725.         (void) flog_delete (plp->pl_rcpz);
  726.         plp->pl_rcpx = plp->pl_rcpy = plp->pl_rcpz = 0;
  727.     }
  728.     npltlines = 0;
  729. }
  730.  
  731. /* let operator select the fields he wants to plot.
  732.  * register them with flog and keep rcfpack() in pltlines[] array.
  733.  */
  734. static
  735. plt_select_fields()
  736. {
  737.     static char hlp[] = "move and RETURN to select a field, or q to quit";
  738.     static char sry[] = "Sorry; can not log any more fields.";
  739.     int r = R_UT, c = C_UTV; /* TODO: start where main was? */
  740.     char buf[64];
  741.     int rcp;
  742.     int i;
  743.  
  744.     plt_reset();
  745.     for (i = 0; i < MAXPLTLINES; i++) {
  746.         sprintf (buf, "select x field for line %d", i+1);
  747.         rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
  748.         if (!rcp)
  749.         break;
  750.         pltlines[i].pl_rcpx = rcp;
  751.         if (flog_add (rcp) < 0) {
  752.         f_msg (sry);
  753.         break;
  754.         }
  755.  
  756.         sprintf (buf, "select y field for line %d", i+1);
  757.         r = unpackr (rcp); c = unpackc (rcp);
  758.         rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
  759.         if (!rcp) {
  760.         (void) flog_delete (pltlines[i].pl_rcpx);
  761.         break;
  762.         }
  763.         if (flog_add (rcp) < 0) {
  764.         (void) flog_delete (pltlines[i].pl_rcpx);
  765.         f_msg (sry);
  766.         break;
  767.         }
  768.         pltlines[i].pl_rcpy = rcp;
  769.  
  770.         sprintf (buf, "select z field for line %d", i+1);
  771.         r = unpackr (rcp); c = unpackc (rcp);
  772.         rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
  773.         if (rcp) {
  774.         if (flog_add (rcp) < 0) {
  775.             (void) flog_delete (pltlines[i].pl_rcpx);
  776.             (void) flog_delete (pltlines[i].pl_rcpy);
  777.             f_msg (sry);
  778.             break;
  779.         }
  780.         pltlines[i].pl_rcpz = rcp;
  781.         r = unpackr (rcp); c = unpackc (rcp);
  782.         }
  783.  
  784.         do {
  785.         sprintf (buf, "enter a one-character label for line %d: ", i+1);
  786.         f_prompt (buf);
  787.         } while (read_line (buf, 1) != 1);
  788.         pltlines[i].pl_label = *buf;
  789.     }
  790.     npltlines = i;
  791. }
  792.  
  793. static
  794. plt_turn_off ()
  795. {
  796.     fclose (plt_fp);
  797.     plt_fp = 0;
  798.     plot_prstate(0);
  799. }
  800.  
  801. /* turn on plotting.
  802.  * establish a file to use (and thereby set plt_fp, the plotting_is_on flag).
  803.  * also check that there is a srch function if it is being plotted.
  804.  */
  805. static
  806. plt_turn_on ()
  807. {
  808.     int sf = rcfpack(R_SRCH, C_SRCH, 0);
  809.     char fn[FNLEN], fnq[64];
  810.     char *optype;
  811.     int n;
  812.     PltLine *plp;
  813.  
  814.     /* insure there is a valid srch function if we are to plot it */
  815.     for (plp = &pltlines[npltlines]; --plp >= pltlines; )
  816.         if ((plp->pl_rcpx == sf || plp->pl_rcpy == sf || plp->pl_rcpz == sf)
  817.             && !prog_isgood()) {
  818.         f_msg ("Plotting search function but it is not defined.");
  819.         return;
  820.         }
  821.  
  822.     /* prompt for file name, giving current as default */
  823.     sprintf (fnq, "file to write <%s>: ", plt_filename);
  824.     f_prompt (fnq);
  825.     n = read_line (fn, sizeof(fn)-1);
  826.  
  827.     /* leave plotting off if type END.
  828.      * reuse same fn if just type \n
  829.      */
  830.     if (n < 0)
  831.         return;
  832.     if (n > 0)
  833.         strcpy (plt_filename, fn);
  834.  
  835.     /* give option to append if file already exists */
  836.     optype = "w";
  837.     if (access (plt_filename, 2) == 0) {
  838.         while (1) {
  839.         f_prompt ("files exists; append or overwrite (a/o)?: ");
  840.         n = read_char();
  841.         if (n == 'a') {
  842.             optype = "a";
  843.             break;
  844.         }
  845.         if (n == 'o')
  846.             break;
  847.         }
  848.     }
  849.  
  850.     /* plotting is on if file opens ok */
  851.     plt_fp = fopen (plt_filename, optype);
  852.     if (!plt_fp) {
  853.         char buf[NC];
  854.         sprintf (buf, "can not open %s: %s", plt_filename, errsys);
  855.         f_prompt (buf);
  856.         (void)read_char();
  857.     }
  858.     plot_prstate (0);
  859. }
  860.  
  861. /* ask operator for a file to plot. if it's ok, do it.
  862.  */
  863. static
  864. plt_file ()
  865. {
  866.     char fn[FNLEN], fnq[64];
  867.     FILE *pfp;
  868.     int n;
  869.  
  870.     /* prompt for file name, giving current as default */
  871.     sprintf (fnq, "file to read <%s>: ", plt_filename);
  872.     f_prompt (fnq);
  873.     n = read_line (fn, sizeof(fn)-1);
  874.  
  875.     /* forget it if type END.
  876.      * reuse same fn if just type \n
  877.      */
  878.     if (n < 0)
  879.         return;
  880.     if (n > 0)
  881.         strcpy (plt_filename, fn);
  882.  
  883.     /* do the plot if file opens ok */
  884.     pfp = fopen (plt_filename, "r");
  885.     if (pfp) {
  886.         if (plt_in_polar)
  887.         plot_polar (pfp);
  888.         else
  889.         plot_cartesian (pfp);
  890.         fclose (pfp);
  891.     } else {
  892.         char buf[NC];
  893.         sprintf (buf, "can not open %s: %s", plt_filename, errsys);
  894.         f_prompt (buf);
  895.         (void)read_char();
  896.     }
  897. }
  898.  
  899. /* plot the given file on the screen in cartesian coords.
  900.  * TODO: add z tags somehow
  901.  * N.B. do whatever you like but redraw the screen when done.
  902.  */
  903. static
  904. plot_cartesian (pfp)
  905. FILE *pfp;
  906. {
  907.     static char fmt[] = "%c,%F,%F";
  908.     double x, y;    /* N.B. be sure these match what your scanf's %F wants*/
  909.     double minx, maxx, miny, maxy;
  910.     char buf[128];
  911.     int npts = 0;
  912.     char c;
  913.  
  914.     /* find ranges and number of points */
  915.     while (fgets (buf, sizeof(buf), pfp)) {
  916.         sscanf (buf, fmt, &c, &x, &y);
  917.         if (npts++ == 0) {
  918.         maxx = minx = x;
  919.         maxy = miny = y;
  920.         } else {
  921.         if (x > maxx) maxx = x;
  922.         else if (x < minx) minx = x;
  923.         if (y > maxy) maxy = y;
  924.         else if (y < miny) miny = y;
  925.         }
  926.     }
  927.  
  928. #define    SMALL    (1e-10)
  929.     if (npts < 2 || fabs(minx-maxx) < SMALL || fabs(miny-maxy) < SMALL)
  930.         f_prompt ("At least two different points required to plot.");
  931.     else {
  932.         /* read file again, this time plotting */
  933.         rewind (pfp);
  934.         c_erase();
  935.         while (fgets (buf, sizeof(buf), pfp)) {
  936.         int row, col;
  937.         sscanf (buf, fmt, &c, &x, &y);
  938.         row = NR-(int)((NR-1)*(y-miny)/(maxy-miny)+0.5);
  939.         col =  1+(int)((NC-1)*(x-minx)/(maxx-minx)+0.5);
  940.         if (row == NR && col == NC)
  941.             col--;    /* avoid lower right scrolling corner */
  942.         f_char (row, col, c);
  943.         }
  944.  
  945.         /* label axes */
  946.         f_double (1, 1, "%.2g", maxy);
  947.         f_double (NR-1, 1, "%.2g", miny);
  948.         f_double (NR, 1, "%.2g", minx);
  949.         f_double (NR, NC-10, "%.2g", maxx);
  950.     }
  951.  
  952.     /* hit any key to resume... */
  953.     (void) read_char();
  954.     redraw_screen (2);    /* full redraw */
  955. }
  956.  
  957. /* plot the given file on the screen in polar coords.
  958.  * first numberic field in plot file is r, second is theta in degrees.
  959.  * TODO: add z tags somehow
  960.  * N.B. do whatever you like but redraw the screen when done.
  961.  */
  962. static
  963. plot_polar (pfp)
  964. FILE *pfp;
  965. {
  966.     static char fmt[] = "%c,%F,%F";
  967.     double r, th;    /* N.B. be sure these match what your scanf's %F wants*/
  968.     double maxr;
  969.     char buf[128];
  970.     int npts = 0;
  971.     char c;
  972.  
  973.     /* find ranges and number of points */
  974.     while (fgets (buf, sizeof(buf), pfp)) {
  975.         sscanf (buf, fmt, &c, &r, &th);
  976.         if (npts++ == 0)
  977.         maxr = r;
  978.         else
  979.         if (r > maxr)
  980.             maxr = r;
  981.     }
  982.  
  983.     if (npts < 2)
  984.         f_prompt ("At least two points required to plot.");
  985.     else {
  986.         /* read file again, this time plotting */
  987.         rewind (pfp);
  988.         c_erase();
  989.         while (fgets (buf, sizeof(buf), pfp)) {
  990.         int row, col;
  991.         double x, y;
  992.         sscanf (buf, fmt, &c, &r, &th);
  993.         x = r * cos(th/57.2958);    /* degs to rads */
  994.         y = r * sin(th/57.2958);
  995.         row = NR-(int)((NR-1)*(y+maxr)/(2.0*maxr)+0.5);
  996.         col =  1+(int)((NC-1)*(x+maxr)/(2.0*maxr)/ASPECT+0.5);
  997.         if (row == NR && col == NC)
  998.             col--;    /* avoid lower right scrolling corner */
  999.         f_char (row, col, c);
  1000.         }
  1001.  
  1002.         /* label radius */
  1003.         f_double (NR/2, NC-10, "%.4g", maxr);
  1004.     }
  1005.  
  1006.     /* hit any key to resume... */
  1007.     (void) read_char();
  1008.     redraw_screen (2);    /* full redraw */
  1009. }
  1010. xXx
  1011. echo x popup.c
  1012. cat > popup.c << 'xXx'
  1013. /* put up a one-line menu consisting of the given fields and let op move
  1014.  * between them with the same methods as sel_fld().
  1015.  * return index of which he picked, or -1 if hit END.
  1016.  */
  1017.  
  1018. #include <stdio.h>
  1019. #include "screen.h"
  1020.  
  1021. #define    FLDGAP    2    /* inter-field gap */
  1022. #define    MAXFLDS    32    /* max number of fields we can handle */
  1023.  
  1024. static char pup[] = "Select: ";
  1025.  
  1026. /* put up an array of strings on prompt line and let op pick one.
  1027.  * start with field fn.
  1028.  * N.B. we do not do much error/bounds checking.
  1029.  */
  1030. popup (fields, fn, nfields)
  1031. char *fields[];
  1032. int fn;
  1033. int nfields;
  1034. {
  1035.     int fcols[MAXFLDS];    /* column to use for each field */
  1036.     int i;
  1037.  
  1038.     if (nfields > MAXFLDS)
  1039.         return (-1);
  1040.  
  1041.     again:
  1042.     /* erase the prompt line; we are going to take it over */
  1043.     c_pos (R_PROMPT, C_PROMPT);
  1044.     c_eol();
  1045.  
  1046.     /* compute starting column for each field */
  1047.     fcols[0] = sizeof(pup);
  1048.     for (i = 1; i < nfields; i++)
  1049.         fcols[i] = fcols[i-1] + strlen (fields[i-1]) + FLDGAP;
  1050.  
  1051.     /* draw each field, with comma after all but last */
  1052.     c_pos (R_PROMPT, 1);
  1053.     fputs (pup, stdout);
  1054.     for (i = 0; i < nfields; i++) {
  1055.         c_pos (R_PROMPT, fcols[i]);
  1056.         printf (i < nfields-1 ? "%s," : "%s", fields[i]);
  1057.     }
  1058.  
  1059.     /* let op choose one now; begin at fn.
  1060.      */
  1061.     i = fn;
  1062.     while (1) {
  1063.         c_pos (R_PROMPT, fcols[i]);
  1064.         switch (read_char()) {
  1065.         case END: return (-1);
  1066.         case REDRAW: redraw_screen(2); goto again;
  1067.         case VERSION: version(); goto again;
  1068.         case '\r': return (i);
  1069.         case 'h':
  1070.         if (--i < 0)
  1071.             i = nfields - 1;
  1072.         break;
  1073.         case 'l':
  1074.         if (++i >= nfields)
  1075.             i = 0;
  1076.         break;
  1077.         }
  1078.     }
  1079. }
  1080. xXx
  1081. echo x precess.c
  1082. cat > precess.c << 'xXx'
  1083. #include <stdio.h>
  1084. #include <math.h>
  1085. #include "astro.h"
  1086.  
  1087. /* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
  1088.  * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
  1089.  * N.B. ra and dec are modifed IN PLACE.
  1090.  */
  1091. precess (mjd1, mjd2, ra, dec)
  1092. double mjd1, mjd2;    /* initial and final epoch modified JDs */
  1093. double *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  1094. {
  1095.     static double lastmjd1, lastmjd2;
  1096.     static double m, n, nyrs;
  1097.     double dra, ddec;    /* ra and dec corrections */
  1098.  
  1099.     if (mjd1 != lastmjd1 || mjd2 != lastmjd2) {
  1100.         double t1, t2; /* Julian centuries of 36525 days since Jan .5 1900 */
  1101.         double m1, n1; /* "constants" for t1 */
  1102.         double m2, n2; /* "constants" for t2 */
  1103.         t1 = mjd1/36525.;
  1104.         t2 = mjd2/36525.;
  1105.         m1 = 3.07234+(1.86e-3*t1);
  1106.         n1 = 20.0468-(8.5e-3*t1);
  1107.         m2 = 3.07234+(1.86e-3*t2);
  1108.         n2 = 20.0468-(8.5e-3*t2);
  1109.         m = (m1+m2)/2;    /* average m for the two epochs */
  1110.         n = (n1+n2)/2;    /* average n for the two epochs */
  1111.         nyrs = (mjd2-mjd1)/365.2425;
  1112.         lastmjd1 = mjd1;
  1113.         lastmjd2 = mjd2;
  1114.     }
  1115.  
  1116.     dra = (m+(n*sin(*ra)*tan(*dec)/15))*7.272205e-5*nyrs;
  1117.     ddec = n*cos(*ra)*4.848137e-6*nyrs;
  1118.     *ra += dra;
  1119.     *dec += ddec;
  1120.     range (ra, 2*PI);
  1121. }
  1122. xXx
  1123. echo x refract.c
  1124. cat > refract.c << 'xXx'
  1125. #include <stdio.h>
  1126. #include <math.h>
  1127. #include "astro.h"
  1128.  
  1129. /* correct the true altitude, ta, for refraction to the apparent altitude, aa,
  1130.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  1131.  * the temperature, tr, in degrees C.
  1132.  */
  1133. refract (pr, tr, ta, aa)
  1134. double pr, tr;
  1135. double ta;
  1136. double *aa;
  1137. {
  1138.     double r;    /* refraction correction*/
  1139.  
  1140.         if (ta >= degrad(15.)) {
  1141.         /* model for altitudes at least 15 degrees above horizon */
  1142.             r = 7.888888e-5*pr/((273+tr)*tan(ta));
  1143.     } else if (ta > degrad(-5.)) {
  1144.         /* hairier model for altitudes at least -5 and below 15 degrees */
  1145.         double a, b, tadeg = raddeg(ta);
  1146.         a = ((2e-5*tadeg+1.96e-2)*tadeg+1.594e-1)*pr;
  1147.         b = (273+tr)*((8.45e-2*tadeg+5.05e-1)*tadeg+1);
  1148.         r = degrad(a/b);
  1149.     } else {
  1150.         /* do nothing if more than 5 degrees below horizon.
  1151.          */
  1152.         r = 0;
  1153.     }
  1154.  
  1155.     *aa  =  ta + r;
  1156. }
  1157.  
  1158. /* correct the apparent altitude, aa, for refraction to the true altitude, ta,
  1159.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  1160.  * the temperature, tr, in degrees C.
  1161.  */
  1162. unrefract (pr, tr, aa, ta)
  1163. double pr, tr;
  1164. double aa;
  1165. double *ta;
  1166. {
  1167.     double err;
  1168.     double appar;
  1169.     double true;
  1170.  
  1171.     /* iterative solution: search for the true that refracts to the
  1172.      *   given apparent.
  1173.      * since refract() is discontinuous at -5 degrees, there is a range
  1174.      *   of apparent altitudes between about -4.5 and -5 degrees that are
  1175.      *   not invertable (the graph of ap vs. true has a vertical step at
  1176.      *   true = -5). thus, the iteration just oscillates if it gets into
  1177.      *   this region. if this happens the iteration is forced to abort.
  1178.      *   of course, this makes unrefract() discontinuous too.
  1179.      */
  1180.     true = aa;
  1181.     do {
  1182.         refract (pr, tr, true, &appar);
  1183.         err = appar - aa;
  1184.         true -= err;
  1185.     } while (fabs(err) >= 1e-6 && true > degrad(-5));
  1186.  
  1187.     *ta = true;
  1188. }
  1189. xXx
  1190. echo x riset.c
  1191. cat > riset.c << 'xXx'
  1192. #include <stdio.h>
  1193. #include <math.h>
  1194. #include "astro.h"
  1195.  
  1196. /* given the true geocentric ra and dec of an object, in radians, the observer's
  1197.  *   latitude in radians, and a horizon displacement correction, dis, in radians
  1198.  *   find the local sidereal times and azimuths of rising and setting, lstr/s
  1199.  *   and azr/s, in radians, respectively.
  1200.  * dis is the vertical displacement from the true position of the horizon. it
  1201.  *   is positive if the apparent position is higher than the true position.
  1202.  *   said another way, it is positive if the shift causes the object to spend
  1203.  *   longer above the horizon. for example, atmospheric refraction is typically
  1204.  *   assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
  1205.  *   would then take on the value +9.89e-3 (radians). On the other hand, if
  1206.  *   your horizon has hills such that your apparent horizon is, say, 1 degree
  1207.  *   above sea level, you would allow for this by setting dis to -1.75e-2
  1208.  *   (radians).
  1209.  * status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
  1210.  */
  1211. riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
  1212. double ra, dec;
  1213. double lat, dis;
  1214. double *lstr, *lsts;
  1215. double *azr, *azs;
  1216. int *status;
  1217. {
  1218.     static double lastlat = 0, slat = 0, clat = 1.0;
  1219.     double dec1, sdec, cdec, tdec;
  1220.     double psi, spsi, cpsi;
  1221.     double h, dh, ch;    /* hr angle, delta and cos */
  1222.  
  1223.     /* avoid needless sin/cos since latitude doesn't change often */
  1224.         if (lat != lastlat) {
  1225.         clat = cos(lat);
  1226.         slat = sin(lat);
  1227.         lastlat = lat;
  1228.     }
  1229.  
  1230.     /* can't cope with objects very near the celestial pole nor if we 
  1231.      * are located very near the earth's poles.
  1232.      */
  1233.     cdec = cos(dec);
  1234.         if (fabs(cdec*clat) < 1e-10) {
  1235.         /* trouble */
  1236.         *status = 2;
  1237.         return;
  1238.     }
  1239.  
  1240.         cpsi = slat/cdec;
  1241.         if (cpsi>1.0) cpsi = 1.0;
  1242.         else if (cpsi<-1.0) cpsi = -1.0;
  1243.         psi = acos(cpsi);
  1244.     spsi = sin(psi);
  1245.  
  1246.         dh = dis*spsi;
  1247.     dec1 = dec + (dis*cpsi);
  1248.         sdec = sin(dec1);
  1249.     tdec = tan(dec1);
  1250.  
  1251.         ch = (-1*slat*tdec)/clat;
  1252.  
  1253.         if (ch < -1) {
  1254.         /* circumpolar */
  1255.         *status = -1;
  1256.         return;
  1257.     }
  1258.         if (ch > 1) {
  1259.         /* never rises */
  1260.         *status = 1;
  1261.         return;
  1262.     }
  1263.  
  1264.         *status = 0;
  1265.     h = acos(ch)+dh;
  1266.  
  1267.         *lstr = 24+radhr(ra-h);
  1268.     *lsts = radhr(ra+h);
  1269.     range (lstr, 24.0);
  1270.     range (lsts, 24.0);
  1271.  
  1272.     *azr = acos(sdec/clat);
  1273.     range (azr, 2*PI);
  1274.         *azs = 2*PI - *azr;
  1275. }
  1276. xXx
  1277. echo x riset_c.c
  1278. cat > riset_c.c << 'xXx'
  1279. /* find rise and set circumstances, ie, riset_cir() and related functions. */
  1280.  
  1281. #include <stdio.h>
  1282. #include <math.h>
  1283. #include "astro.h"
  1284. #include "circum.h"
  1285. #include "screen.h"    /* just for SUN and MOON */
  1286.  
  1287. #define    TRACE(x)    {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
  1288.  
  1289. #define    STDREF    degrad(34./60.)    /* nominal horizon refraction amount */
  1290. #define    TWIREF    degrad(18.)    /* twilight horizon displacement */
  1291.  
  1292. /* find where and when a body, p, will rise and set and
  1293.  *   its transit circumstances. all times are local, angles rads e of n.
  1294.  * return 0 if just returned same stuff as previous call, else 1 if new.
  1295.  * status is set from the RS_* #defines in circum.h.
  1296.  * also used to find astro twilight by calling with dis of 18 degrees.
  1297.  */
  1298. riset_cir (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
  1299. int p;        /* one of the body defines in astro.h or screen.h */
  1300. Now *np;
  1301. int hzn;    /* STDHZN or ADPHZN or TWILIGHT */
  1302. double *ltr, *lts; /* local rise and set times */
  1303. double *ltt;    /* local transit time */
  1304. double *azr, *azs; /* local rise and set azimuths, rads e of n */
  1305. double *altt;    /* local altitude at transit */
  1306. int *status;    /* one or more of the RS_* defines */
  1307. {
  1308.     typedef struct {
  1309.         Now l_now;
  1310.         double l_ltr, l_lts, l_ltt, l_azr, l_azs, l_altt;
  1311.         int l_hzn;
  1312.         int l_status;
  1313.     } Last;
  1314.     /* must be in same order as the astro.h/screen.h #define's */
  1315.     static Last last[NOBJ];
  1316.     Last *lp;
  1317.     int new;
  1318.  
  1319.     lp = last + p;
  1320.     if (same_cir (np, &lp->l_now) && same_lday (np, &lp->l_now)
  1321.                         && lp->l_hzn == hzn) {
  1322.         *ltr = lp->l_ltr;
  1323.         *lts = lp->l_lts;
  1324.         *ltt = lp->l_ltt;
  1325.         *azr = lp->l_azr;
  1326.         *azs = lp->l_azs;
  1327.         *altt = lp->l_altt;
  1328.         *status = lp->l_status;
  1329.         new = 0;
  1330.     } else {
  1331.         *status = 0;
  1332.         iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status);
  1333.         lp->l_ltr = *ltr;
  1334.         lp->l_lts = *lts;
  1335.         lp->l_ltt = *ltt;
  1336.         lp->l_azr = *azr;
  1337.         lp->l_azs = *azs;
  1338.         lp->l_altt = *altt;
  1339.         lp->l_status = *status;
  1340.         lp->l_hzn = hzn;
  1341.         lp->l_now = *np;
  1342.         new = 1;
  1343.     }
  1344.     return (new);
  1345. }
  1346.  
  1347. static
  1348. iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
  1349. int p;
  1350. Now *np;
  1351. int hzn;
  1352. double *ltr, *lts, *ltt;    /* local times of rise, set and transit */
  1353. double *azr, *azs, *altt;/* local azimuths of rise, set and transit altitude */
  1354. int *status;
  1355. {
  1356. #define    MAXPASSES    6
  1357.     double lstr, lsts, lstt;    /* local sidereal times of rising/setting */
  1358.     double mjd0;        /* mjd estimates of rise/set event */
  1359.     double lnoon;        /* mjd of local noon */
  1360.     double x;        /* discarded tmp value */
  1361.     Now n;            /* just used to call now_lst() */
  1362.     double lst;        /* lst at local noon */
  1363.     double diff, lastdiff;    /* iterative improvement to mjd0 */
  1364.     int pass;
  1365.     int rss;
  1366.  
  1367.     /* first approximation is to find rise/set times of a fixed object
  1368.      * in its position at local noon.
  1369.      */
  1370.     lnoon = mjd_day(mjd - tz/24.0) + (12.0 + tz)/24.0; /*mjd of local noon*/
  1371.     n.n_mjd = lnoon;
  1372.     n.n_lng = lng;
  1373.     now_lst (&n, &lst);    /* lst at local noon */
  1374.     mjd0 = lnoon;
  1375.     stationary_riset (p,mjd0,np,hzn,&lstr,&lsts,&lstt,&x,&x,&x,&rss);
  1376.     chkrss:
  1377.     switch (rss) {
  1378.     case  0:  break;
  1379.     case  1: *status = RS_NEVERUP; return;
  1380.     case -1: *status = RS_CIRCUMPOLAR; goto transit;
  1381.     default: *status = RS_ERROR; return;
  1382.     }
  1383.  
  1384.     /* find a better approximation to the rising circumstances based on
  1385.      * more passes, each using a "fixed" object at the location at
  1386.      * previous approximation of the rise time.
  1387.      */
  1388.     lastdiff = 1000.0;
  1389.     for (pass = 1; pass < MAXPASSES; pass++) {
  1390.         diff = (lstr - lst)*SIDRATE; /* next guess at rise time wrt noon */
  1391.         if (diff > 12.0)
  1392.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  1393.         else if (diff < -12.0)
  1394.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  1395.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of rise */
  1396.         stationary_riset (p,mjd0,np,hzn,&lstr,&x,&x,azr,&x,&x,&rss);
  1397.         if (rss != 0) goto chkrss;
  1398.         if (fabs (diff - lastdiff) < 0.25/60.0)
  1399.         break;             /* converged to better than 15 secs */
  1400.         lastdiff = diff;
  1401.     }
  1402.     if (pass == MAXPASSES)
  1403.         *status |= RS_NORISE;    /* didn't converge - no rise today */
  1404.     else {
  1405.         *ltr = 12.0 + diff;
  1406.         if (p != MOON &&
  1407.             (*ltr <= 24.0*(1.0-SIDRATE) || *ltr >= 24.0*SIDRATE))
  1408.         *status |= RS_2RISES;
  1409.     }
  1410.  
  1411.     /* find a better approximation to the setting circumstances based on
  1412.      * more passes, each using a "fixed" object at the location at
  1413.      * previous approximation of the set time.
  1414.      */
  1415.     lastdiff = 1000.0;
  1416.     for (pass = 1; pass < MAXPASSES; pass++) {
  1417.         diff = (lsts - lst)*SIDRATE; /* next guess at set time wrt noon */
  1418.         if (diff > 12.0)
  1419.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  1420.         else if (diff < -12.0)
  1421.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  1422.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of set */
  1423.         stationary_riset (p,mjd0,np,hzn,&x,&lsts,&x,&x,azs,&x,&rss);
  1424.         if (rss != 0) goto chkrss;
  1425.         if (fabs (diff - lastdiff) < 0.25/60.0)
  1426.         break;             /* converged to better than 15 secs */
  1427.         lastdiff = diff;
  1428.     }
  1429.     if (pass == MAXPASSES)
  1430.         *status |= RS_NOSET;    /* didn't converge - no set today */
  1431.     else {
  1432.         *lts = 12.0 + diff;
  1433.         if (p != MOON &&
  1434.             (*lts <= 24.0*(1.0-SIDRATE) || *lts >= 24.0*SIDRATE))
  1435.         *status |= RS_2SETS;
  1436.     }
  1437.  
  1438.     transit:
  1439.     /* find a better approximation to the transit circumstances based on
  1440.      * more passes, each using a "fixed" object at the location at
  1441.      * previous approximation of the transit time.
  1442.      */
  1443.     lastdiff = 1000.0;
  1444.     for (pass = 1; pass < MAXPASSES; pass++) {
  1445.         diff = (lstt - lst)*SIDRATE; /*next guess at transit time wrt noon*/
  1446.         if (diff > 12.0)
  1447.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  1448.         else if (diff < -12.0)
  1449.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  1450.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of set */
  1451.         stationary_riset (p,mjd0,np,hzn,&x,&x,&lstt,&x,&x,altt,&rss);
  1452.         if (fabs (diff - lastdiff) < 0.25/60.0)
  1453.         break;             /* converged to better than 15 secs */
  1454.         lastdiff = diff;
  1455.     }
  1456.     if (pass == MAXPASSES || rss != 0)
  1457.         *status |= RS_NOTRANS;    /* didn't converge - no transit today */
  1458.     else {
  1459.         *ltt = 12.0 + diff;
  1460.         if (p != MOON &&
  1461.             (*ltt <= 24.0*(1.0-SIDRATE) || *ltt >= 24.0*SIDRATE))
  1462.         *status |= RS_2TRANS;
  1463.     }
  1464. }
  1465.  
  1466. static
  1467. stationary_riset (p, mjd0, np, hzn, lstr, lsts, lstt, azr, azs, altt, status)
  1468. int p;
  1469. double mjd0;
  1470. Now *np;
  1471. int hzn;
  1472. double *lstr, *lsts, *lstt;
  1473. double *azr, *azs, *altt;
  1474. int *status;
  1475. {
  1476.     double dis;
  1477.     Now n;
  1478.     Sky s;
  1479.  
  1480.     switch (hzn) {
  1481.     case STDHZN:
  1482.         /* nominal atmospheric refraction.
  1483.          * then add nominal moon or sun semi-diameter, as appropriate.
  1484.          */
  1485.         dis = STDREF;
  1486.         break;
  1487.     case TWILIGHT:
  1488.         if (p != SUN) {
  1489.         c_erase();
  1490.         printf ("BUG 1!\n\r");
  1491.         exit(1);
  1492.         }
  1493.         dis = TWIREF;
  1494.         break;
  1495.     default:
  1496.         /* horizon displacement is refraction plus object's semi-diameter */
  1497.         unrefract (pressure, temp, 0.0, &dis);
  1498.         dis = -dis;
  1499.         n = *np;
  1500.         n.n_mjd = mjd0;
  1501.         (void) body_cir (p, 0.0, &n, &s);
  1502.         dis += degrad(s.s_size/3600./2.0);
  1503.     }
  1504.  
  1505.     switch (p) {
  1506.     case SUN:
  1507.         if (hzn == STDHZN)
  1508.         dis += degrad (32./60./2.);
  1509.         fixedsunriset (mjd0,np,dis,lstr,lsts,lstt,azr,azs,altt,status);
  1510.         break;
  1511.     case MOON:
  1512.         if (hzn == STDHZN)
  1513.         dis += degrad (32./60./2.);
  1514.         fixedmoonriset (mjd0,np,dis,lstr,lsts,lstt,azr,azs,altt,status);
  1515.         break;
  1516.     default:
  1517.         if (hzn == STDHZN) {
  1518.         /* assume planets have negligible diameters.
  1519.          * also, need to set s if hzn was STDHZN.
  1520.          */
  1521.         n = *np;
  1522.         n.n_mjd = mjd0;
  1523.         (void) body_cir (p, 0.0, &n, &s);
  1524.         }
  1525.         riset (s.s_ra, s.s_dec, lat, dis, lstr, lsts, azr, azs, status);
  1526.         transit (s.s_ra, s.s_dec, np, lstt, altt);
  1527.     }
  1528. }
  1529.  
  1530. /* find the local rise/set sidereal times and azimuths of an object fixed
  1531.  * at the ra/dec of the sun on the given mjd time as seen from lat.
  1532.  * times are for the upper limb. dis should account for refraction and
  1533.  *   sun's semi-diameter; we add corrections for nutation and aberration.
  1534.  */
  1535. static
  1536. fixedsunriset (mjd0, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
  1537. double mjd0;
  1538. Now *np;
  1539. double dis;
  1540. double *lstr, *lsts, *lstt;
  1541. double *azr, *azs, *altt;
  1542. int *status;
  1543. {
  1544.     double lsn, rsn;
  1545.     double deps, dpsi;
  1546.     double r, d;
  1547.  
  1548.     /* find ecliptic position of sun at mjd */
  1549.     sunpos (mjd0, &lsn, &rsn);
  1550.  
  1551.     /* allow for 20.4" aberation and nutation */
  1552.     nutation (mjd0, &deps, &dpsi);
  1553.         lsn += dpsi - degrad(20.4/3600.0);
  1554.  
  1555.     /* convert ecliptic to equatorial coords */
  1556.     ecl_eq (mjd0, 0.0, lsn, &r, &d);
  1557.  
  1558.     /* find circumstances for given horizon displacement */
  1559.     riset (r, d, lat, dis, lstr, lsts, azr, azs, status);
  1560.     transit (r, d, np, lstt, altt);
  1561. }
  1562.  
  1563. /* find the local sidereal rise/set times and azimuths of an object fixed
  1564.  * at the ra/dec of the moon on the given mjd time as seen from lat.
  1565.  * times are for the upper limb. dis should account for refraction and moon's
  1566.  *    semi-diameter; we add corrections for parallax, and nutation.
  1567.  * accuracy is to nearest minute.
  1568.  */
  1569. static
  1570. fixedmoonriset (mjd0, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
  1571. double mjd0;
  1572. Now *np;
  1573. double dis;
  1574. double *lstr, *lsts, *lstt;
  1575. double *azr, *azs, *altt;
  1576. int *status;
  1577. {
  1578.     double lam, bet, hp;
  1579.     double deps, dpsi;
  1580.     double r, d;
  1581.     double ha;
  1582.  
  1583.     /* find geocentric ecliptic location and equatorial horizontal parallax
  1584.      * of moon at mjd.
  1585.      */
  1586.     moon (mjd0, &lam, &bet, &hp);
  1587.  
  1588.     /* allow for nutation */
  1589.     nutation (mjd0, &deps, &dpsi);
  1590.         lam += dpsi;
  1591.  
  1592.     /* convert ecliptic to equatorial coords */
  1593.     ecl_eq (mjd0, bet, lam, &r, &d);
  1594.  
  1595.     /* find local sidereal times of rise/set times/azimuths for given
  1596.      * equatorial coords, allowing for
  1597.      * .27249*sin(hp)   parallax (hp is radius of earth from moon;
  1598.      *                  .27249 is radius of moon from earth where
  1599.      *                  the ratio is the dia_moon/dia_earth).
  1600.      * hp               nominal angular earth radius. subtract because
  1601.      *                  tangential line-of-sight makes moon appear lower
  1602.      *                  as move out from center of earth.
  1603.      */
  1604.     dis += .27249*sin(hp) - hp;
  1605.     riset (r, d, lat, dis, lstr, lsts, azr, azs, status);
  1606.  
  1607.     /* account for parallax on the meridian (0 hour-angle).
  1608.      * TODO: is this really right? other stuff too? better way?? help!
  1609.      */
  1610.     ta_par (0.0, d, lat, height, hp, &ha, &d);
  1611.     transit (r+ha, d, np, lstt, altt);
  1612. }
  1613.  
  1614. /* find when and how hi object at (r,d) is when it transits. */
  1615. static
  1616. transit (r, d, np, lstt, altt)
  1617. double r, d;    /* ra and dec, rads */
  1618. Now *np;    /* for refraction info */
  1619. double *lstt;    /* local sidereal time of transit */
  1620. double *altt;    /* local, refracted, altitude at time of transit */
  1621. {
  1622.     *lstt = radhr(r);
  1623.     *altt = PI/2 - lat + d;
  1624.     if (*altt > PI/2)
  1625.         *altt = PI - *altt;
  1626.     refract (pressure, temp, *altt, altt);
  1627. }
  1628. xXx
  1629.