This is precess.c in view mode; [Download] [Up]
#include <stdio.h>
#include <math.h>
#include "astro.h"
#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)))
/* 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.
*/
/*
* 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
*/
precess (mjd1, mjd2, ra, dec)
double mjd1, mjd2; /* initial and final epoch modified JDs */
double *ra, *dec; /* ra/dec for mjd1 in, for mjd2 out */
{
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;
mjd_year (mjd1, &from_equinox);
mjd_year (mjd2, &to_equinox);
alpha_in = raddeg(*ra);
delta_in = raddeg(*dec);
/* From from_equinox to 2000.0 */
if (from_equinox != 2000.0) {
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 (to_equinox != 2000.0) {
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);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.