diff --git a/libastro/astro.h b/libastro/astro.h index c7b4967..6557fad 100644 --- a/libastro/astro.h +++ b/libastro/astro.h @@ -4,12 +4,18 @@ #include #ifndef PI -#define PI 3.141592653589793 +#define PI 3.141592653589793238 +#endif + +#ifndef TWOPI +#define TWOPI 6.283185307179586477 #endif /* conversions among hours (of ra), degrees and radians. */ #define degrad(x) ((x)*PI/180.) #define raddeg(x) ((x)*180./PI) +#define radarcsec(x) ((x)*2.06264806247096355e5) +#define arcsecrad(x) ((x)*4.84813681109535993e-6) #define hrdeg(x) ((x)*15.) #define deghr(x) ((x)/15.) #define hrrad(x) degrad(hrdeg(x)) diff --git a/libastro/precess.c b/libastro/precess.c index 2816708..dd0847a 100644 --- a/libastro/precess.c +++ b/libastro/precess.c @@ -1,144 +1,69 @@ #include #include - #include "astro.h" -static void precess_hiprec (double mjd1, double mjd2, double *ra, double *dec); +/* Precession computation to or from J2000 + * If dir < 0 go from equinox to J2000 + * If dir > 0 go from J2000 to equinox + */ +static void +do_precess ( +double equinox, double *alpha, double *delta, int dir +) +{ + double T; + double zeta_A, z_A, theta_A; + double A, B, C; + /* Precession progresses about 1 arc second in .047 years = 17 days + * Threshold must be at most 1.7 days for 0.1 arcsec accuracy */ + double threshhold=1; -#define DCOS(x) cos(degrad(x)) -#define DSIN(x) sin(degrad(x)) -#define DASIN(x) raddeg(asin(x)) -#define DATAN2(y,x) raddeg(atan2((y),(x))) + if (fabs (equinox - J2000) < threshhold) + return; -/* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2. - * the epochs are given by their modified JDs, mjd1 and mjd2, respectively. - * N.B. ra and dec are modifed IN PLACE. + T = (equinox - J2000)/36525.0; + + /* Compute coefficients in arcseconds from the Naval + * Observatory, Her Majesty's Nautical Amanac Office; + * Rutherford Appleton Laboratory, ISSN 0737-6421. */ + + zeta_A = arcsecrad(2.650545 + T*(2306.083227+T*(0.2988499 + T*(0.01801828 +T*(-5.971e-6-3.173e-7*T))))); + z_A = arcsecrad(-2.650545+T*(2306.077181+T*(1.0927348+T*(0.01826837 +T*(-28.596e-6-2.904e-7*T))))); + theta_A = arcsecrad(T*(2004.191903-T*(0.4294934+T*(0.04182264+T*(7.089e-6+1.274e-7*T))))); + + if (dir<0){ + /* If dir is negative, go from equinox to 2000.0 */ + A = sin(*alpha - z_A) * cos(*delta); + B = cos(*alpha - z_A) * cos(theta_A) * cos(*delta) + + sin(theta_A) * sin(*delta); + C = -cos(*alpha - z_A) * sin(theta_A) * cos(*delta) + + cos(theta_A) * sin(*delta); + *alpha = atan2(A,B) - zeta_A; + } else { + /* If dir is positive, go from 2000.0 to equinox */ + A = sin(*alpha + zeta_A) * cos(*delta); + B = cos(*alpha + zeta_A) * cos(theta_A) * cos(*delta) + - sin(theta_A) * sin(*delta); + C = cos(*alpha + zeta_A) * sin(theta_A) * cos(*delta) + + cos(theta_A) * sin(*delta); + *alpha = atan2(A,B) + z_A; + } + range (alpha, TWOPI); + *delta = asin(C); +} + +/* Corrects ra and dec, both in radians, for precession from epoch 1 to + * epoch 2. The epochs are given by their modified JDs, from_mjd and + * to_mjd, respectively. N.B. ra and dec are modifed IN PLACE. */ void precess ( -double mjd1, double mjd2, /* initial and final epoch modified JDs */ -double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */ +double from_mjd, double to_mjd, /* initial and final epoch modified JDs */ +double *ra, double *dec) /* ra/dec for from_mjd in, for to_mjd out */ { - precess_hiprec (mjd1, mjd2, ra, dec); + /* From mjd1 to J2000, "backwards" inplace transformation : */ + do_precess(from_mjd, ra, dec, -1); + /* From J2000 to mjd2, "forward" inplace transformation : */ + do_precess(to_mjd, ra, dec, 1); } - -/* - * Copyright (c) 1990 by Craig Counterman. All rights reserved. - * - * This software may be redistributed freely, not sold. - * This copyright notice and disclaimer of warranty must remain - * unchanged. - * - * No representation is made about the suitability of this - * software for any purpose. It is provided "as is" without express or - * implied warranty, to the extent permitted by applicable law. - * - * Rigorous precession. From Astronomical Ephemeris 1989, p. B18 - * - * 96-06-20 Hayo Hase : theta_a corrected - */ -static void -precess_hiprec ( -double mjd1, double mjd2, /* initial and final epoch modified JDs */ -double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */ -{ - static double last_mjd1 = -213.432, last_from; - static double last_mjd2 = -213.432, last_to; - double zeta_A, z_A, theta_A; - double T; - double A, B, C; - double alpha, delta; - double alpha_in, delta_in; - double from_equinox, to_equinox; - double alpha2000, delta2000; - - /* convert mjds to years; - * avoid the remarkably expensive calls to mjd_year() - */ - if (last_mjd1 == mjd1) - from_equinox = last_from; - else { - mjd_year (mjd1, &from_equinox); - last_mjd1 = mjd1; - last_from = from_equinox; - } - if (last_mjd2 == mjd2) - to_equinox = last_to; - else { - mjd_year (mjd2, &to_equinox); - last_mjd2 = mjd2; - last_to = to_equinox; - } - - /* convert coords in rads to degs */ - alpha_in = raddeg(*ra); - delta_in = raddeg(*dec); - - /* precession progresses about 1 arc second in .047 years */ - /* From from_equinox to 2000.0 */ - if (fabs (from_equinox-2000.0) > .02) { - T = (from_equinox - 2000.0)/100.0; - zeta_A = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T; - z_A = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T; - theta_A = 0.5567530* T - 0.0001185* T*T - 0.0000116* T*T*T; - - A = DSIN(alpha_in - z_A) * DCOS(delta_in); - B = DCOS(alpha_in - z_A) * DCOS(theta_A) * DCOS(delta_in) - + DSIN(theta_A) * DSIN(delta_in); - C = -DCOS(alpha_in - z_A) * DSIN(theta_A) * DCOS(delta_in) - + DCOS(theta_A) * DSIN(delta_in); - - alpha2000 = DATAN2(A,B) - zeta_A; - range (&alpha2000, 360.0); - delta2000 = DASIN(C); - } else { - /* should get the same answer, but this could improve accruacy */ - alpha2000 = alpha_in; - delta2000 = delta_in; - }; - - - /* From 2000.0 to to_equinox */ - if (fabs (to_equinox - 2000.0) > .02) { - T = (to_equinox - 2000.0)/100.0; - zeta_A = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T; - z_A = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T; - theta_A = 0.5567530* T - 0.0001185* T*T - 0.0000116* T*T*T; - - A = DSIN(alpha2000 + zeta_A) * DCOS(delta2000); - B = DCOS(alpha2000 + zeta_A) * DCOS(theta_A) * DCOS(delta2000) - - DSIN(theta_A) * DSIN(delta2000); - C = DCOS(alpha2000 + zeta_A) * DSIN(theta_A) * DCOS(delta2000) - + DCOS(theta_A) * DSIN(delta2000); - - alpha = DATAN2(A,B) + z_A; - range(&alpha, 360.0); - delta = DASIN(C); - } else { - /* should get the same answer, but this could improve accruacy */ - alpha = alpha2000; - delta = delta2000; - }; - - *ra = degrad(alpha); - *dec = degrad(delta); -} - -#if 0 -static void -precess_fast ( -double mjd1, double mjd2, /* initial and final epoch modified JDs */ -double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */ -{ -#define N degrad (20.0468/3600.0) -#define M hrrad (3.07234/3600.0) - double nyrs; - - nyrs = (mjd2 - mjd1)/365.2425; - *dec += N * cos(*ra) * nyrs; - *ra += (M + (N * sin(*ra) * tan(*dec))) * nyrs; - range (ra, 2.0*PI); -} -#endif -