home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part19 / plans.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  18.6 KB  |  626 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5.  
  6. #if defined(__STDC__) || defined(__cplusplus)
  7. #define P_(s) s
  8. #else
  9. #define P_(s) ()
  10. #endif
  11.  
  12. extern void anomaly P_((double ma, double s, double *nu, double *ea));
  13. extern void pelement P_((double Mjd, double plan[8][9]));
  14. extern void range P_((double *v, double r));
  15. extern void sunpos P_((double Mjd, double *lsn, double *rsn));
  16.  
  17. void plans P_((double mjd, int p, double *lpd0, double *psi0, double *rp0, double *rho0, double *lam, double *bet, double *dia, double *mag));
  18. static void aux_jsun P_((double t, double *x1, double *x2, double *x3, double *x4, double *x5, double *x6));
  19. static void masun P_((double mjd, double *mas));
  20. static void p_mercury P_((double map[], double *dl, double *dr));
  21. static void p_venus P_((double t, double mas, double map[], double *dl, double *dr, double *dml, double *dm));
  22. static void p_mars P_((double mas, double map[], double *dl, double *dr, double *dml, double *dm));
  23. static void p_jupiter P_((double t, double s, double *dml, double *ds, double *dm, double *da));
  24. static void p_saturn P_((double t, double s, double *dml, double *ds, double *dm, double *da, double *dhl));
  25. static void p_uranus P_((double t, double s, double *dl, double *dr, double *dml, double *ds, double *dm, double *da, double *dhl));
  26. static void p_neptune P_((double t, double s, double *dl, double *dr, double *dml, double *ds, double *dm, double *da, double *dhl));
  27.  
  28. #undef P_
  29.  
  30. #define    TWOPI        (2*PI)
  31. #define    mod2PI(x)    ((x) - (long)((x)/TWOPI)*TWOPI)
  32.  
  33. /* given a modified Julian date, mjd, and a planet, p, find:
  34.  *   lpd0: heliocentric longitude, 
  35.  *   psi0: heliocentric latitude,
  36.  *   rp0:  distance from the sun to the planet, 
  37.  *   rho0: distance from the Earth to the planet,
  38.  *         none corrected for light time, ie, they are the true values for the
  39.  *         given instant.
  40.  *   lam:  geocentric ecliptic longitude, 
  41.  *   bet:  geocentric ecliptic latitude,
  42.  *         each corrected for light time, ie, they are the apparent values as
  43.  *       seen from the center of the Earth for the given instant.
  44.  *   dia:  angular diameter in arcsec at 1 AU, 
  45.  *   mag:  visual magnitude when 1 AU from sun and earth at 0 phase angle.
  46.  *
  47.  * all angles are in radians, all distances in AU.
  48.  * the mean orbital elements are found by calling pelement(), then mutual
  49.  *   perturbation corrections are applied as necessary.
  50.  *
  51.  * corrections for nutation and abberation must be made by the caller. The RA 
  52.  *   and DEC calculated from the fully-corrected ecliptic coordinates are then
  53.  *   the apparent geocentric coordinates. Further corrections can be made, if
  54.  *   required, for atmospheric refraction and geocentric parallax although the
  55.  *   intrinsic error herein of about 10 arcseconds is usually the dominant
  56.  *   error at this stage.
  57.  * TODO: combine the several intermediate expressions when get a good compiler.
  58.  */
  59. void
  60. plans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
  61. double mjd;
  62. int p;
  63. double *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
  64. {
  65.     static double plan[8][9];
  66.     static double lastmjd = -10000;
  67.     double dl;    /* perturbation correction for longitude */
  68.     double dr;    /*  "   orbital radius */
  69.     double dml;    /*  "   mean longitude */
  70.     double ds;    /*  "   eccentricity */
  71.     double dm;    /*  "   mean anomaly */
  72.     double da;    /*  "   semi-major axis */
  73.     double dhl;    /*  "   heliocentric longitude */
  74.     double lsn, rsn;/* true geocentric longitude of sun and sun-earth rad */
  75.     double mas;    /* mean anomaly of the sun */
  76.     double re;    /* radius of earth's orbit */
  77.     double lg;    /* longitude of earth */
  78.     double map[8];    /* array of mean anomalies for each planet */
  79.     double lpd, psi, rp, rho;
  80.     double ll, sll, cll;
  81.     double t;
  82.     double dt;
  83.     int pass;
  84.     int j;
  85.     double s, ma;
  86.     double nu, ea;
  87.     double lp, om;
  88.     double lo, slo, clo;
  89.     double inc, y;
  90.     double spsi, cpsi;
  91.     double rpd;
  92.  
  93.     /* only need to fill in plan[] once for a given mjd */
  94.     if (mjd != lastmjd) {
  95.         pelement (mjd, plan);
  96.         lastmjd = mjd;
  97.     }
  98.  
  99.     dt = 0;
  100.     t = mjd/36525.;
  101.     sunpos (mjd, &lsn, &rsn);
  102.     masun (mjd, &mas);
  103.         re = rsn;
  104.     lg = lsn+PI;
  105.  
  106.     /* first find the true position of the planet at mjd.
  107.      * then repeat a second time for a slightly different time based
  108.      * on the position found in the first pass to account for light-travel
  109.      * time.
  110.      */
  111.     for (pass = 0; pass < 2; pass++) {
  112.  
  113.         for (j = 0; j < 8; j++)
  114.         map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]);
  115.  
  116.         /* set initial corrections to 0.
  117.          * then modify as necessary for the planet of interest.
  118.          */
  119.         dl = 0;
  120.         dr = 0;
  121.         dml = 0;
  122.         ds = 0;
  123.         dm = 0;
  124.         da = 0;
  125.         dhl = 0;
  126.  
  127.         switch (p) {
  128.  
  129.         case MERCURY:
  130.         p_mercury (map, &dl, &dr);
  131.         break;
  132.  
  133.         case VENUS:
  134.         p_venus (t, mas, map, &dl, &dr, &dml, &dm);
  135.         break;
  136.  
  137.         case MARS:
  138.         p_mars (mas, map, &dl, &dr, &dml, &dm);
  139.         break;
  140.  
  141.         case JUPITER:
  142.         p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da);
  143.         break;
  144.  
  145.         case SATURN:
  146.         p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl);
  147.         break;
  148.  
  149.         case URANUS:
  150.         p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
  151.         break;
  152.  
  153.         case NEPTUNE:
  154.         p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
  155.         break;
  156.  
  157.         case PLUTO:
  158.         /* no perturbation theory for pluto */
  159.         break;
  160.         }
  161.  
  162.         s = plan[p][3]+ds;
  163.         ma = map[p]+dm;
  164.         anomaly (ma, s, &nu, &ea);
  165.         rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu));
  166.         lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm);
  167.         lp = degrad(lp);
  168.         om = degrad(plan[p][5]);
  169.         lo = lp-om;
  170.         slo = sin(lo);
  171.         clo = cos(lo);
  172.         inc = degrad(plan[p][4]);
  173.         rp = rp+dr;
  174.         spsi = slo*sin(inc);
  175.         y = slo*cos(inc);
  176.         psi = asin(spsi)+dhl;
  177.         spsi = sin(psi);
  178.         lpd = atan(y/clo)+om+degrad(dl);
  179.         if (clo<0) lpd += PI;
  180.         range (&lpd, TWOPI);
  181.         cpsi = cos(psi);
  182.         rpd = rp*cpsi;
  183.         ll = lpd-lg;
  184.         rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll));
  185.  
  186.         /* when we view a planet we see it in the position it occupied
  187.          * dt days ago, where rho is the distance between it and earth,
  188.          * in AU. use this as the new time for the next pass.
  189.          */
  190.         dt = rho*5.775518e-3;
  191.  
  192.         if (pass == 0) {
  193.         /* save heliocentric coordinates after first pass since, being
  194.          * true, they are NOT to be corrected for light-travel time.
  195.          */
  196.         *lpd0 = lpd;
  197.         range (lpd0, TWOPI);
  198.         *psi0 = psi;
  199.         *rp0 = rp;
  200.         *rho0 = rho;
  201.         }
  202.     }
  203.  
  204.         sll = sin(ll);
  205.     cll = cos(ll);
  206.         if (p < MARS) 
  207.         *lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI;
  208.     else
  209.         *lam = atan(re*sll/(rpd-re*cll))+lpd;
  210.     range (lam, TWOPI);
  211.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll));
  212.     *dia = plan[p][7];
  213.     *mag = plan[p][8];
  214. }
  215.  
  216. /* set auxilliary variables used for jupiter, saturn, uranus, and neptune */
  217. static void
  218. aux_jsun (t, x1, x2, x3, x4, x5, x6)
  219. double t;
  220. double *x1, *x2, *x3, *x4, *x5, *x6;
  221. {
  222.         *x1 = t/5+0.1;
  223.         *x2 = mod2PI(4.14473+5.29691e1*t);
  224.         *x3 = mod2PI(4.641118+2.132991e1*t);
  225.         *x4 = mod2PI(4.250177+7.478172*t);
  226.         *x5 = 5 * *x3 - 2 * *x2;
  227.     *x6 = 2 * *x2 - 6 * *x3 + 3 * *x4;
  228. }
  229.  
  230. /* find the mean anomaly of the sun at mjd.
  231.  * this is the same as that used in sun() but when it was converted to C it
  232.  * was not known it would be required outside that routine.
  233.  * TODO: add an argument to sun() to return mas and eliminate this routine.
  234.  */
  235. static void
  236. masun (mjd, mas)
  237. double mjd;
  238. double *mas;
  239. {
  240.     double t, t2;
  241.     double a, b;
  242.  
  243.     t = mjd/36525;
  244.     t2 = t*t;
  245.     a = 9.999736042e1*t;
  246.     b = 360.*(a-(long)a);
  247.     *mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b);
  248. }
  249.  
  250. /* perturbations for mercury */
  251. static void
  252. p_mercury (map, dl, dr)
  253. double map[];
  254. double *dl, *dr;
  255. {
  256.     *dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+
  257.          1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+
  258.          9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+
  259.          7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1);
  260.  
  261.     *dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+
  262.          6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+
  263.          5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+
  264.          3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699);
  265. }
  266.  
  267. /* ....venus */
  268. static void
  269. p_venus (t, mas, map, dl, dr, dml, dm)
  270. double t, mas, map[];
  271. double *dl, *dr, *dml, *dm;
  272. {
  273.     *dml = degrad (7.7e-4*sin(4.1406+t*2.6227));
  274.     *dm = *dml;
  275.  
  276.     *dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+
  277.          1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+
  278.          1.36e-3*cos(mas-map[2-1]-2.0788)+
  279.          9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+
  280.          8.2e-4*cos(map[3]-map[2-1]-3.6318);
  281.  
  282.     *dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+
  283.          1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+
  284.          6.887e-6*cos(map[3]-map[2-1]-2.06106)+
  285.          5.172e-6*cos(mas-map[2-1]-5.08065e-1)+
  286.          3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+
  287.          3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+
  288.          3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1);
  289. }
  290.  
  291. /* ....mars */
  292. static void
  293. p_mars (mas, map, dl, dr, dml, dm)
  294. double mas, map[];
  295. double *dl, *dr, *dml, *dm;
  296. {
  297.     double a;
  298.  
  299.     a = 3*map[3]-8*map[2]+4*mas;
  300.     *dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
  301.     *dm = *dml;
  302.  
  303.     *dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
  304.          6.07e-3*cos(2*map[3]-map[2]-3.2873)+
  305.          4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
  306.          3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
  307.          2.38e-3*cos(mas-map[2]+6.1256e-1)+
  308.          2.04e-3*cos(2*mas-3*map[2]+2.7688)+
  309.          1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
  310.          1.36e-3*cos(2*mas-4*map[2]+2.6894)+
  311.          1.04e-3*cos(map[3]+3.0749e-1);
  312.  
  313.     *dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
  314.          5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
  315.          3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
  316.          1.5996e-5*cos(mas-map[2]-9.69618e-1)+
  317.          1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
  318.          8.966e-6*cos(map[3]-2*map[2]+7.61225e-1);
  319.      *dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
  320.          7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
  321.          6.62e-6*cos(mas-2*map[2]+1.97575)+
  322.          4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
  323.          4.693e-6*cos(3*mas-5*map[2]+3.32665)+
  324.          4.571e-6*cos(2*mas-4*map[2]+4.27086)+
  325.          4.409e-6*cos(3*map[3]-map[2]-2.02158);
  326. }
  327.  
  328. /* ....jupiter */
  329. static void
  330. p_jupiter (t, s, dml, ds, dm, da)
  331. double t, s;
  332. double *dml, *ds, *dm, *da;
  333. {
  334.     double dp;
  335.     double x1, x2, x3, x4, x5, x6, x7;
  336.     double sx3, cx3, s2x3, c2x3;
  337.         double sx5, cx5, s2x5;
  338.     double sx6;
  339.         double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7;
  340.  
  341.     aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
  342.         x7 = x3-x2;
  343.     sx3 = sin(x3);
  344.     cx3 = cos(x3);
  345.         s2x3 = sin(2*x3);
  346.     c2x3 = cos(2*x3);
  347.         sx5 = sin(x5);
  348.     cx5 = cos(x5);
  349.         s2x5 = sin(2*x5);
  350.     sx6 = sin(x6);
  351.         sx7 = sin(x7);
  352.     cx7 = cos(x7);
  353.         s2x7 = sin(2*x7);
  354.     c2x7 = cos(2*x7);
  355.         s3x7 = sin(3*x7);
  356.     c3x7 = cos(3*x7);
  357.         s4x7 = sin(4*x7);
  358.     c4x7 = cos(4*x7);
  359.         c5x7 = cos(5*x7);
  360.  
  361.     *dml = (3.31364e-1-(1.0281e-2+4.692e-3*x1)*x1)*sx5+
  362.           (3.228e-3-(6.4436e-2-2.075e-3*x1)*x1)*cx5-
  363.           (3.083e-3+(2.75e-4-4.89e-4*x1)*x1)*s2x5+
  364.           2.472e-3*sx6+1.3619e-2*sx7+1.8472e-2*s2x7+6.717e-3*s3x7+
  365.           2.775e-3*s4x7+6.417e-3*s2x7*sx3+
  366.           (7.275e-3-1.253e-3*x1)*sx7*sx3+
  367.           2.439e-3*s3x7*sx3-(3.5681e-2+1.208e-3*x1)*sx7*cx3;
  368.         *dml += -3.767e-3*c2x7*sx3-(3.3839e-2+1.125e-3*x1)*cx7*sx3-
  369.           4.261e-3*s2x7*cx3+
  370.           (1.161e-3*x1-6.333e-3)*cx7*cx3+
  371.           2.178e-3*cx3-6.675e-3*c2x7*cx3-2.664e-3*c3x7*cx3-
  372.           2.572e-3*sx7*s2x3-3.567e-3*s2x7*s2x3+2.094e-3*cx7*c2x3+
  373.           3.342e-3*c2x7*c2x3;
  374.     *dml = degrad(*dml);
  375.  
  376.     *ds = (3606+(130-43*x1)*x1)*sx5+(1289-580*x1)*cx5-6764*sx7*sx3-
  377.          1110*s2x7*sx3-224*s3x7*sx3-204*sx3+(1284+116*x1)*cx7*sx3+
  378.          188*c2x7*sx3+(1460+130*x1)*sx7*cx3+224*s2x7*cx3-817*cx3+
  379.          6074*cx3*cx7+992*c2x7*cx3+
  380.          508*c3x7*cx3+230*c4x7*cx3+108*c5x7*cx3;
  381.     *ds += -(956+73*x1)*sx7*s2x3+448*s2x7*s2x3+137*s3x7*s2x3+
  382.          (108*x1-997)*cx7*s2x3+480*c2x7*s2x3+148*c3x7*s2x3+
  383.          (99*x1-956)*sx7*c2x3+490*s2x7*c2x3+
  384.          158*s3x7*c2x3+179*c2x3+(1024+75*x1)*cx7*c2x3-
  385.          437*c2x7*c2x3-132*c3x7*c2x3;
  386.     *ds *= 1e-7;
  387.  
  388.     dp = (7.192e-3-3.147e-3*x1)*sx5-4.344e-3*sx3+
  389.          (x1*(1.97e-4*x1-6.75e-4)-2.0428e-2)*cx5+
  390.          3.4036e-2*cx7*sx3+(7.269e-3+6.72e-4*x1)*sx7*sx3+
  391.          5.614e-3*c2x7*sx3+2.964e-3*c3x7*sx3+3.7761e-2*sx7*cx3+
  392.          6.158e-3*s2x7*cx3-
  393.          6.603e-3*cx7*cx3-5.356e-3*sx7*s2x3+2.722e-3*s2x7*s2x3+
  394.          4.483e-3*cx7*s2x3-2.642e-3*c2x7*s2x3+4.403e-3*sx7*c2x3-
  395.          2.536e-3*s2x7*c2x3+5.547e-3*cx7*c2x3-2.689e-3*c2x7*c2x3;
  396.  
  397.     *dm = *dml-(degrad(dp)/s);
  398.  
  399.     *da = 205*cx7-263*cx5+693*c2x7+312*c3x7+147*c4x7+299*sx7*sx3+
  400.          181*c2x7*sx3+204*s2x7*cx3+111*s3x7*cx3-337*cx7*cx3-
  401.          111*c2x7*cx3;
  402.     *da *= 1e-6;
  403. }
  404.  
  405. /* ....saturn */
  406. static void
  407. p_saturn (t, s, dml, ds, dm, da, dhl)
  408. double t, s;
  409. double *dml, *ds, *dm, *da, *dhl;
  410. {
  411.     double dp;
  412.     double x1, x2, x3, x4, x5, x6, x7, x8;
  413.     double sx3, cx3, s2x3, c2x3, s3x3, c3x3, s4x3, c4x3;
  414.         double sx5, cx5, s2x5, c2x5;
  415.     double sx6;
  416.         double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7, s5x7;
  417.     double s2x8, c2x8, s3x8, c3x8;
  418.  
  419.     aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
  420.         x7 = x3-x2;
  421.     sx3 = sin(x3);
  422.     cx3 = cos(x3);
  423.         s2x3 = sin(2*x3);
  424.     c2x3 = cos(2*x3);
  425.         sx5 = sin(x5);
  426.     cx5 = cos(x5);
  427.         s2x5 = sin(2*x5);
  428.     sx6 = sin(x6);
  429.         sx7 = sin(x7);
  430.     cx7 = cos(x7);
  431.         s2x7 = sin(2*x7);
  432.     c2x7 = cos(2*x7);
  433.         s3x7 = sin(3*x7);
  434.     c3x7 = cos(3*x7);
  435.         s4x7 = sin(4*x7);
  436.     c4x7 = cos(4*x7);
  437.         c5x7 = cos(5*x7);
  438.  
  439.     s3x3 = sin(3*x3);
  440.     c3x3 = cos(3*x3);
  441.     s4x3 = sin(4*x3);
  442.     c4x3 = cos(4*x3);
  443.     c2x5 = cos(2*x5);
  444.     s5x7 = sin(5*x7);
  445.     x8 = x4-x3;
  446.     s2x8 = sin(2*x8);
  447.     c2x8 = cos(2*x8);
  448.     s3x8 = sin(3*x8);
  449.     c3x8 = cos(3*x8);
  450.  
  451.     *dml = 7.581e-3*s2x5-7.986e-3*sx6-1.48811e-1*sx7-4.0786e-2*s2x7-
  452.           (8.14181e-1-(1.815e-2-1.6714e-2*x1)*x1)*sx5-
  453.           (1.0497e-2-(1.60906e-1-4.1e-3*x1)*x1)*cx5-1.5208e-2*s3x7-
  454.           6.339e-3*s4x7-6.244e-3*sx3-1.65e-2*s2x7*sx3+
  455.           (8.931e-3+2.728e-3*x1)*sx7*sx3-5.775e-3*s3x7*sx3+
  456.           (8.1344e-2+3.206e-3*x1)*cx7*sx3+1.5019e-2*c2x7*sx3;
  457.     *dml += (8.5581e-2+2.494e-3*x1)*sx7*cx3+1.4394e-2*c2x7*cx3+
  458.           (2.5328e-2-3.117e-3*x1)*cx7*cx3+
  459.           6.319e-3*c3x7*cx3+6.369e-3*sx7*s2x3+9.156e-3*s2x7*s2x3+
  460.           7.525e-3*s3x8*s2x3-5.236e-3*cx7*c2x3-7.736e-3*c2x7*c2x3-
  461.           7.528e-3*c3x8*c2x3;
  462.     *dml = degrad(*dml);
  463.  
  464.     *ds = (-7927+(2548+91*x1)*x1)*sx5+(13381+(1226-253*x1)*x1)*cx5+
  465.          (248-121*x1)*s2x5-(305+91*x1)*c2x5+412*s2x7+12415*sx3+
  466.          (390-617*x1)*sx7*sx3+(165-204*x1)*s2x7*sx3+26599*cx7*sx3-
  467.          4687*c2x7*sx3-1870*c3x7*sx3-821*c4x7*sx3-
  468.          377*c5x7*sx3+497*c2x8*sx3+(163-611*x1)*cx3;
  469.     *ds += -12696*sx7*cx3-4200*s2x7*cx3-1503*s3x7*cx3-619*s4x7*cx3-
  470.          268*s5x7*cx3-(282+1306*x1)*cx7*cx3+(-86+230*x1)*c2x7*cx3+
  471.          461*s2x8*cx3-350*s2x3+(2211-286*x1)*sx7*s2x3-
  472.          2208*s2x7*s2x3-568*s3x7*s2x3-346*s4x7*s2x3-
  473.          (2780+222*x1)*cx7*s2x3+(2022+263*x1)*c2x7*s2x3+248*c3x7*s2x3+
  474.          242*s3x8*s2x3+467*c3x8*s2x3-490*c2x3-(2842+279*x1)*sx7*c2x3;
  475.     *ds += (128+226*x1)*s2x7*c2x3+224*s3x7*c2x3+
  476.          (-1594+282*x1)*cx7*c2x3+(2162-207*x1)*c2x7*c2x3+
  477.          561*c3x7*c2x3+343*c4x7*c2x3+469*s3x8*c2x3-242*c3x8*c2x3-
  478.          205*sx7*s3x3+262*s3x7*s3x3+208*cx7*c3x3-271*c3x7*c3x3-
  479.          382*c3x7*s4x3-376*s3x7*c4x3;
  480.     *ds *= 1e-7;
  481.  
  482.     dp = (7.7108e-2+(7.186e-3-1.533e-3*x1)*x1)*sx5-7.075e-3*sx7+
  483.          (4.5803e-2-(1.4766e-2+5.36e-4*x1)*x1)*cx5-7.2586e-2*cx3-
  484.          7.5825e-2*sx7*sx3-2.4839e-2*s2x7*sx3-8.631e-3*s3x7*sx3-
  485.          1.50383e-1*cx7*cx3+2.6897e-2*c2x7*cx3+1.0053e-2*c3x7*cx3-
  486.          (1.3597e-2+1.719e-3*x1)*sx7*s2x3+1.1981e-2*s2x7*c2x3;
  487.     dp += -(7.742e-3-1.517e-3*x1)*cx7*s2x3+
  488.          (1.3586e-2-1.375e-3*x1)*c2x7*c2x3-
  489.          (1.3667e-2-1.239e-3*x1)*sx7*c2x3+
  490.          (1.4861e-2+1.136e-3*x1)*cx7*c2x3-
  491.          (1.3064e-2+1.628e-3*x1)*c2x7*c2x3;
  492.  
  493.     *dm = *dml-(degrad(dp)/s);
  494.  
  495.     *da = 572*sx5-1590*s2x7*cx3+2933*cx5-647*s3x7*cx3+33629*cx7-
  496.          344*s4x7*cx3-3081*c2x7+2885*cx7*cx3-1423*c3x7+
  497.          (2172+102*x1)*c2x7*cx3-671*c4x7+296*c3x7*cx3-320*c5x7-
  498.          267*s2x7*s2x3+1098*sx3-778*cx7*s2x3-2812*sx7*sx3;
  499.     *da += 495*c2x7*s2x3+688*s2x7*sx3+250*c3x7*s2x3-393*s3x7*sx3-
  500.          856*sx7*c2x3-228*s4x7*sx3+441*s2x7*c2x3+2138*cx7*sx3+
  501.          296*c2x7*c2x3-999*c2x7*sx3+211*c3x7*c2x3-642*c3x7*sx3-
  502.          427*sx7*s3x3-325*c4x7*sx3+398*s3x7*s3x3-890*cx3+
  503.          344*cx7*c3x3+2206*sx7*cx3-427*c3x7*c3x3;
  504.     *da *= 1e-6;
  505.  
  506.     *dhl = 7.47e-4*cx7*sx3+1.069e-3*cx7*cx3+2.108e-3*s2x7*s2x3+
  507.           1.261e-3*c2x7*s2x3+1.236e-3*s2x7*c2x3-2.075e-3*c2x7*c2x3;
  508.     *dhl = degrad(*dhl);
  509. }
  510.  
  511. /* ....uranus */
  512. static void
  513. p_uranus (t, s, dl, dr, dml, ds, dm, da, dhl)
  514. double t, s;
  515. double *dl, *dr, *dml, *ds, *dm, *da, *dhl;
  516. {
  517.     double dp;
  518.     double x1, x2, x3, x4, x5, x6;
  519.     double x8, x9, x10, x11, x12;
  520.     double sx4, cx4, s2x4, c2x4;
  521.     double sx9, cx9, s2x9, c2x9;
  522.     double sx11, cx11;
  523.  
  524.     aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
  525.  
  526.         x8 = mod2PI(1.46205+3.81337*t);
  527.         x9 = 2*x8-x4;
  528.     sx9 = sin(x9);
  529.     cx9 = cos(x9);
  530.         s2x9 = sin(2*x9);
  531.     c2x9 = cos(2*x9);
  532.  
  533.     x10 = x4-x2;
  534.     x11 = x4-x3;
  535.     x12 = x8-x4;
  536.  
  537.     *dml = (8.64319e-1-1.583e-3*x1)*sx9+(8.2222e-2-6.833e-3*x1)*cx9+
  538.           3.6017e-2*s2x9-3.019e-3*c2x9+8.122e-3*sin(x6);
  539.     *dml = degrad(*dml);
  540.  
  541.     dp = 1.20303e-1*sx9+6.197e-3*s2x9+(1.9472e-2-9.47e-4*x1)*cx9;
  542.     *dm = *dml-(degrad(dp)/s);
  543.  
  544.     *ds = (163*x1-3349)*sx9+20981*cx9+1311*c2x9;
  545.     *ds *= 1e-7;
  546.  
  547.     *da = -3.825e-3*cx9;
  548.  
  549.     *dl = (1.0122e-2-9.88e-4*x1)*sin(x4+x11)+
  550.          (-3.8581e-2+(2.031e-3-1.91e-3*x1)*x1)*cos(x4+x11)+
  551.          (3.4964e-2-(1.038e-3-8.68e-4*x1)*x1)*cos(2*x4+x11)+
  552.          5.594e-3*sin(x4+3*x12)-1.4808e-2*sin(x10)-
  553.          5.794e-3*sin(x11)+2.347e-3*cos(x11)+9.872e-3*sin(x12)+
  554.          8.803e-3*sin(2*x12)-4.308e-3*sin(3*x12);
  555.  
  556.     sx11 = sin(x11);
  557.     cx11 = cos(x11);
  558.     sx4 = sin(x4);
  559.     cx4 = cos(x4);
  560.     s2x4 = sin(2*x4);
  561.     c2x4 = cos(2*x4);
  562.     *dhl = (4.58e-4*sx11-6.42e-4*cx11-5.17e-4*cos(4*x12))*sx4-
  563.           (3.47e-4*sx11+8.53e-4*cx11+5.17e-4*sin(4*x11))*cx4+
  564.           4.03e-4*(cos(2*x12)*s2x4+sin(2*x12)*c2x4);
  565.     *dhl = degrad(*dhl);
  566.  
  567.     *dr = -25948+4985*cos(x10)-1230*cx4+3354*cos(x11)+904*cos(2*x12)+
  568.          894*(cos(x12)-cos(3*x12))+(5795*cx4-1165*sx4+1388*c2x4)*sx11+
  569.          (1351*cx4+5702*sx4+1388*s2x4)*cos(x11);
  570.     *dr *= 1e-6;
  571. }
  572.  
  573. /* ....neptune */
  574. static void
  575. p_neptune (t, s, dl, dr, dml, ds, dm, da, dhl)
  576. double t, s;
  577. double *dl, *dr, *dml, *ds, *dm, *da, *dhl;
  578. {
  579.     double dp;
  580.     double x1, x2, x3, x4, x5, x6;
  581.     double x8, x9, x10, x11, x12;
  582.     double sx8, cx8;
  583.     double sx9, cx9, s2x9, c2x9;
  584.     double s2x12, c2x12;
  585.  
  586.     aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
  587.  
  588.         x8 = mod2PI(1.46205+3.81337*t);
  589.         x9 = 2*x8-x4;
  590.     sx9 = sin(x9);
  591.     cx9 = cos(x9);
  592.         s2x9 = sin(2*x9);
  593.     c2x9 = cos(2*x9);
  594.  
  595.     x10 = x8-x2;
  596.     x11 = x8-x3;
  597.     x12 = x8-x4;
  598.  
  599.     *dml = (1.089e-3*x1-5.89833e-1)*sx9+(4.658e-3*x1-5.6094e-2)*cx9-
  600.           2.4286e-2*s2x9;
  601.     *dml = degrad(*dml);
  602.  
  603.     dp = 2.4039e-2*sx9-2.5303e-2*cx9+6.206e-3*s2x9-5.992e-3*c2x9;
  604.  
  605.     *dm = *dml-(degrad(dp)/s);
  606.  
  607.     *ds = 4389*sx9+1129*s2x9+4262*cx9+1089*c2x9;
  608.     *ds *= 1e-7;
  609.  
  610.     *da = 8189*cx9-817*sx9+781*c2x9;
  611.     *da *= 1e-6;
  612.  
  613.     s2x12 = sin(2*x12);
  614.     c2x12 = cos(2*x12);
  615.     sx8 = sin(x8);
  616.     cx8 = cos(x8);
  617.     *dl = -9.556e-3*sin(x10)-5.178e-3*sin(x11)+2.572e-3*s2x12-
  618.          2.972e-3*c2x12*sx8-2.833e-3*s2x12*cx8;
  619.  
  620.     *dhl = 3.36e-4*c2x12*sx8+3.64e-4*s2x12*cx8;
  621.     *dhl = degrad(*dhl);
  622.  
  623.     *dr = -40596+4992*cos(x10)+2744*cos(x11)+2044*cos(x12)+1051*c2x12;
  624.     *dr *= 1e-6;
  625. }
  626.