mirror of https://github.com/XEphem/XEphem.git
Precession formulae from 2000 Astro. Almanac
Closes #45 by merging in its changes; thank you, Francois Meyer!
This commit is contained in:
francois.meyer
Brandon Rhodes
parent
266ef93744
commit
5522e1039b
|
|
@ -4,12 +4,18 @@
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
|
|
||||||
#ifndef PI
|
#ifndef PI
|
||||||
#define PI 3.141592653589793
|
#define PI 3.141592653589793238
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifndef TWOPI
|
||||||
|
#define TWOPI 6.283185307179586477
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
/* conversions among hours (of ra), degrees and radians. */
|
/* conversions among hours (of ra), degrees and radians. */
|
||||||
#define degrad(x) ((x)*PI/180.)
|
#define degrad(x) ((x)*PI/180.)
|
||||||
#define raddeg(x) ((x)*180./PI)
|
#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 hrdeg(x) ((x)*15.)
|
||||||
#define deghr(x) ((x)/15.)
|
#define deghr(x) ((x)/15.)
|
||||||
#define hrrad(x) degrad(hrdeg(x))
|
#define hrrad(x) degrad(hrdeg(x))
|
||||||
|
|
|
||||||
|
|
@ -1,144 +1,69 @@
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <math.h>
|
#include <math.h>
|
||||||
|
|
||||||
#include "astro.h"
|
#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))
|
if (fabs (equinox - J2000) < threshhold)
|
||||||
#define DSIN(x) sin(degrad(x))
|
return;
|
||||||
#define DASIN(x) raddeg(asin(x))
|
|
||||||
#define DATAN2(y,x) raddeg(atan2((y),(x)))
|
|
||||||
|
|
||||||
/* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
|
T = (equinox - J2000)/36525.0;
|
||||||
* the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
|
|
||||||
* N.B. ra and dec are modifed IN PLACE.
|
/* 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
|
void
|
||||||
precess (
|
precess (
|
||||||
double mjd1, double mjd2, /* initial and final epoch modified JDs */
|
double from_mjd, double to_mjd, /* initial and final epoch modified JDs */
|
||||||
double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */
|
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 <hase@wettzell.ifag.de>: 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
|
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue