home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / reduce.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  1.9 KB  |  78 lines
  1. #include <math.h>
  2. #include "astro.h"
  3.  
  4.  
  5. #if defined(__STDC__) || defined(__cplusplus)
  6. #define P_(s) s
  7. #else
  8. #define P_(s) ()
  9. #endif
  10.  
  11. extern void range P_((double *v, double r));
  12.  
  13. void reduce_elements P_((double mjd0, double mjd, double inc0, double ap0, double om0, double *inc, double *ap, double *om));
  14.  
  15. #undef P_
  16.  
  17. /* convert those orbital elements that change from epoch mjd0 to epoch mjd.
  18.  */
  19. void
  20. reduce_elements (mjd0, mjd, inc0, ap0, om0, inc, ap, om)
  21. double mjd0;    /* initial epoch */
  22. double mjd;    /* desired epoch */
  23. double inc0;    /* initial inclination, rads */
  24. double ap0;    /* initial argument of perihelion, as an mjd */
  25. double om0;    /* initial long of ascending node, rads */
  26. double *inc;    /* desired inclination, rads */
  27. double *ap;    /* desired epoch of perihelion, as an mjd */
  28. double *om;    /* desired long of ascending node, rads */
  29. {
  30.     double t0, t1;
  31.     double tt, tt2, t02, tt3;
  32.     double eta, th, th0;
  33.     double a, b;
  34.     double dap;
  35.     double cinc, sinc;
  36.     double ot, sot, cot, ot1;
  37.     double seta, ceta;
  38.  
  39.     t0 = mjd0/365250.0;
  40.     t1 = mjd/365250.0;
  41.  
  42.     tt = t1-t0;
  43.     tt2 = tt*tt;
  44.         t02 = t0*t0;
  45.     tt3 = tt*tt2;
  46.         eta = (471.07-6.75*t0+.57*t02)*tt+(.57*t0-3.37)*tt2+.05*tt3;
  47.         th0 = 32869.0*t0+56*t02-(8694+55*t0)*tt+3*tt2;
  48.         eta = degrad(eta/3600.0);
  49.         th0 = degrad((th0/3600.0)+173.950833);
  50.         th = (50256.41+222.29*t0+.26*t02)*tt+(111.15+.26*t0)*tt2+.1*tt3;
  51.         th = th0+degrad(th/3600.0);
  52.     cinc = cos(inc0);
  53.         sinc = sin(inc0);
  54.     ot = om0-th0;
  55.     sot = sin(ot);
  56.         cot = cos(ot);
  57.     seta = sin(eta);
  58.         ceta = cos(eta);
  59.     a = sinc*sot;
  60.         b = ceta*sinc*cot-seta*cinc;
  61.     ot1 = atan(a/b);
  62.         if (b<0) ot1 += PI;
  63.         b = sinc*ceta-cinc*seta*cot;
  64.         a = -1*seta*sot;
  65.     dap = atan(a/b);
  66.         if (b<0) dap += PI;
  67.  
  68.         *ap = ap0+dap;
  69.     range (ap, 2*PI);
  70.         *om = ot1+th;
  71.     range (om, 2*PI);
  72.  
  73.         if (inc0<.175)
  74.         *inc = asin(a/sin(dap));
  75.     else
  76.         *inc = 1.570796327-asin((cinc*ceta)+(sinc*seta*cot));
  77. }
  78.