home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part18 / precess.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  4.3 KB  |  151 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4. #include "preferences.h"
  5.  
  6.  
  7. #if defined(__STDC__) || defined(__cplusplus)
  8. #define P_(s) s
  9. #else
  10. #define P_(s) ()
  11. #endif
  12.  
  13. extern void mjd_year P_((double Mjd, double *yr));
  14. extern void range P_((double *v, double r));
  15.  
  16. void precess P_((double mjd1, double mjd2, double *ra, double *dec));
  17. static void precess_hiprec P_((double mjd1, double mjd2, double *ra, double *dec));
  18. static void precess_fast P_((double mjd1, double mjd2, double *ra, double *dec));
  19.  
  20. #undef P_
  21.  
  22. #define    DCOS(x)        cos(degrad(x))
  23. #define    DSIN(x)        sin(degrad(x))
  24. #define    DASIN(x)    raddeg(asin(x))
  25. #define    DATAN2(y,x)    raddeg(atan2((y),(x)))
  26.  
  27. /* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
  28.  * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
  29.  * N.B. ra and dec are modifed IN PLACE.
  30.  */
  31. void
  32. precess (mjd1, mjd2, ra, dec)
  33. double mjd1, mjd2;    /* initial and final epoch modified JDs */
  34. double *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  35. {
  36.     if (pref_get(PREF_ALGO) == PREF_ACCURATE)
  37.         precess_hiprec (mjd1, mjd2, ra, dec);
  38.     else
  39.         precess_fast (mjd1, mjd2, ra, dec);
  40. }
  41.  
  42. /*
  43.  * Copyright (c) 1990 by Craig Counterman. All rights reserved.
  44.  *
  45.  * This software may be redistributed freely, not sold.
  46.  * This copyright notice and disclaimer of warranty must remain
  47.  *    unchanged. 
  48.  *
  49.  * No representation is made about the suitability of this
  50.  * software for any purpose.  It is provided "as is" without express or
  51.  * implied warranty, to the extent permitted by applicable law.
  52.  *
  53.  * Rigorous precession. From Astronomical Ephemeris 1989, p. B18
  54.  */
  55. static void
  56. precess_hiprec (mjd1, mjd2, ra, dec)
  57. double mjd1, mjd2;    /* initial and final epoch modified JDs */
  58. double *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  59. {
  60.     double zeta_A, z_A, theta_A;
  61.     double T;
  62.     double A, B, C;
  63.     double alpha, delta;
  64.     double alpha_in, delta_in;
  65.     double from_equinox, to_equinox;
  66.     double alpha2000, delta2000;
  67.  
  68.     mjd_year (mjd1, &from_equinox);
  69.     mjd_year (mjd2, &to_equinox);
  70.     alpha_in = raddeg(*ra);
  71.     delta_in = raddeg(*dec);
  72.  
  73.     /* From from_equinox to 2000.0 */
  74.     if (from_equinox != 2000.0) {
  75.         T = (from_equinox - 2000.0)/100.0;
  76.         zeta_A  = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
  77.         z_A     = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
  78.         theta_A = 0.5567530* T - 0.0001185* T*T + 0.0000116* T*T*T;
  79.  
  80.         A = DSIN(alpha_in - z_A) * DCOS(delta_in);
  81.         B = DCOS(alpha_in - z_A) * DCOS(theta_A) * DCOS(delta_in)
  82.           + DSIN(theta_A) * DSIN(delta_in);
  83.         C = -DCOS(alpha_in - z_A) * DSIN(theta_A) * DCOS(delta_in)
  84.           + DCOS(theta_A) * DSIN(delta_in);
  85.  
  86.         alpha2000 = DATAN2(A,B) - zeta_A;
  87.         range (&alpha2000, 360.0);
  88.         delta2000 = DASIN(C);
  89.     } else {
  90.         /* should get the same answer, but this could improve accruacy */
  91.         alpha2000 = alpha_in;
  92.         delta2000 = delta_in;
  93.     };
  94.  
  95.  
  96.     /* From 2000.0 to to_equinox */
  97.     if (to_equinox != 2000.0) {
  98.         T = (to_equinox - 2000.0)/100.0;
  99.         zeta_A  = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
  100.         z_A     = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
  101.         theta_A = 0.5567530* T - 0.0001185* T*T + 0.0000116* T*T*T;
  102.  
  103.         A = DSIN(alpha2000 + zeta_A) * DCOS(delta2000);
  104.         B = DCOS(alpha2000 + zeta_A) * DCOS(theta_A) * DCOS(delta2000)
  105.           - DSIN(theta_A) * DSIN(delta2000);
  106.         C = DCOS(alpha2000 + zeta_A) * DSIN(theta_A) * DCOS(delta2000)
  107.           + DCOS(theta_A) * DSIN(delta2000);
  108.  
  109.         alpha = DATAN2(A,B) + z_A;
  110.         range(&alpha, 360.0);
  111.         delta = DASIN(C);
  112.     } else {
  113.         /* should get the same answer, but this could improve accruacy */
  114.         alpha = alpha2000;
  115.         delta = delta2000;
  116.     };
  117.  
  118.     *ra = degrad(alpha);
  119.     *dec = degrad(delta);
  120. }
  121.  
  122. static void
  123. precess_fast (mjd1, mjd2, ra, dec)
  124. double mjd1, mjd2;    /* initial and final epoch modified JDs */
  125. double *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  126. {
  127.     double nyrs;
  128.     double t1, t2;
  129.     double na, nb;
  130.     double ma, mb;
  131.     double n, m;
  132.     double ddec, dra;
  133.  
  134.     nyrs = (mjd2 - mjd1)/365.2425;
  135.     t2 = mjd2/36525.0;
  136.     t1 = mjd1/36525.0;
  137.     nb = 20.0468 - (8.5e-3 * t2);
  138.     na = 20.0468 - (8.5e-3 * t1);
  139.     mb = 3.07234 + (1.86e-3 * t2);
  140.     ma = 3.07234 + (1.86e-3 * t1);
  141.     n = (na + nb)/2.0;
  142.     m = (ma + mb)/2.0;
  143.  
  144.     ddec = n*cos(*ra) * 4.848137e-6 * nyrs;
  145.     dra  = (m + (n * sin(*ra) * tan(*dec)/15.0)) * 7.272205e-5 * nyrs;
  146.  
  147.     *dec += ddec;
  148.     *ra += dra;
  149.     range (ra, 2.0*PI);
  150. }
  151.