home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / ephem421.zip / REDUCE.C < prev    next >
C/C++ Source or Header  |  1990-09-13  |  2KB  |  64 lines
  1. #include <math.h>
  2. #include "astro.h"
  3.  
  4. /* convert those orbital elements that change from epoch mjd0 to epoch mjd.
  5.  */
  6. reduce_elements (mjd0, mjd, inc0, ap0, om0, inc, ap, om)
  7. double mjd0;    /* initial epoch */
  8. double mjd;    /* desired epoch */
  9. double inc0;    /* initial inclination, rads */
  10. double ap0;    /* initial argument of perihelion, as an mjd */
  11. double om0;    /* initial long of ascending node, rads */
  12. double *inc;    /* desired inclination, rads */
  13. double *ap;    /* desired epoch of perihelion, as an mjd */
  14. double *om;    /* desired long of ascending node, rads */
  15. {
  16.     double t0, t1;
  17.     double tt, tt2, t02, tt3;
  18.     double eta, th, th0;
  19.     double a, b;
  20.     double dap;
  21.     double cinc, sinc;
  22.     double ot, sot, cot, ot1;
  23.     double seta, ceta;
  24.  
  25.     t0 = mjd0/365250.0;
  26.     t1 = mjd/365250.0;
  27.  
  28.     tt = t1-t0;
  29.     tt2 = tt*tt;
  30.         t02 = t0*t0;
  31.     tt3 = tt*tt2;
  32.         eta = (471.07-6.75*t0+.57*t02)*tt+(.57*t0-3.37)*tt2+.05*tt3;
  33.         th0 = 32869.0*t0+56*t02-(8694+55*t0)*tt+3*tt2;
  34.         eta = degrad(eta/3600.0);
  35.         th0 = degrad((th0/3600.0)+173.950833);
  36.         th = (50256.41+222.29*t0+.26*t02)*tt+(111.15+.26*t0)*tt2+.1*tt3;
  37.         th = th0+degrad(th/3600.0);
  38.     cinc = cos(inc0);
  39.         sinc = sin(inc0);
  40.     ot = om0-th0;
  41.     sot = sin(ot);
  42.         cot = cos(ot);
  43.     seta = sin(eta);
  44.         ceta = cos(eta);
  45.     a = sinc*sot;
  46.         b = ceta*sinc*cot-seta*cinc;
  47.     ot1 = atan(a/b);
  48.         if (b<0) ot1 += PI;
  49.         b = sinc*ceta-cinc*seta*cot;
  50.         a = -1*seta*sot;
  51.     dap = atan(a/b);
  52.         if (b<0) dap += PI;
  53.  
  54.         *ap = ap0+dap;
  55.     range (ap, 2*PI);
  56.         *om = ot1+th;
  57.     range (om, 2*PI);
  58.  
  59.         if (inc0<.175)
  60.         *inc = asin(a/sin(dap));
  61.     else
  62.         *inc = 1.570796327-asin((cinc*ceta)+(sinc*seta*cot));
  63. }
  64.