home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / MISC / EDUCATIO / EPHEM421.ZIP / PRECESS.C < prev    next >
Encoding:
C/C++ Source or Header  |  1990-09-13  |  1.5 KB  |  49 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
  6.  * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
  7.  * N.B. ra and dec are modifed IN PLACE.
  8.  * TODO: find a better algorithm; this one is not even symmetric.
  9.  */
  10. precess (mjd1, mjd2, ra, dec)
  11. double mjd1, mjd2;    /* initial and final epoch modified JDs */
  12. double *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  13. {
  14.     static double lastmjd1 = -10000, lastmjd2 = -10000;
  15.     static double m, n, nyrs;
  16.     double dra, ddec;    /* ra and dec corrections */
  17.  
  18.     if (mjd1 != lastmjd1 || mjd2 != lastmjd2) {
  19.         double t1, t2; /* Julian centuries of 36525 days since Jan .5 1900*/
  20.         double m1, n1; /* "constants" for t1 */
  21.         double m2, n2; /* "constants" for t2 */
  22.         t1 = mjd1/36525.;
  23.         t2 = mjd2/36525.;
  24.         m1 = 3.07234+(1.86e-3*t1);
  25.         n1 = 20.0468-(8.5e-3*t1);
  26.         m2 = 3.07234+(1.86e-3*t2);
  27.         n2 = 20.0468-(8.5e-3*t2);
  28.         m = (m1+m2)/2;    /* average m for the two epochs */
  29.         n = (n1+n2)/2;    /* average n for the two epochs */
  30.         nyrs = (mjd2-mjd1)/365.2425;
  31.         lastmjd1 = mjd1;
  32.         lastmjd2 = mjd2;
  33.     }
  34.  
  35.     dra = (m+(n*sin(*ra)*tan(*dec)/15))*7.272205e-5*nyrs;
  36.     ddec = n*cos(*ra)*4.848137e-6*nyrs;
  37.     *ra += dra;
  38.     *dec += ddec;
  39.     /* added by ECD */
  40.     if (*dec > PI/2) {
  41.         *dec = PI - *dec;
  42.         *ra += PI;
  43.     } else if (*dec < -PI/2) {
  44.         *dec = -PI - *dec;
  45.         *ra += PI;
  46.     }
  47.     range (ra, 2*PI);
  48. }
  49.