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