home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / EPHEM421.ZIP / 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.689