This is precess.c in view mode; [Download] [Up]
/************************************************/
/* precess.c */
/* */
/* Precession of the equinox and ecliptic */
/* from epoch Julian date J to or from J2000.0 */
/* */
/* Program by Steve Moshier. */
/************************************************/
/***** description
*
* $Id: precess.c,v 1.3 1993/04/21 21:43:02 craig Exp $
*
* IAU Coefficients are from:
* J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
* "Expressions for the Precession Quantities Based upon the IAU
* (1976) System of Astronomical Constants," Astronomy and
* Astrophysics 58, 1-16 (1977).
*
* Newer formulas that cover a much longer time span are from:
* J. Laskar, "Secular terms of classical planetary theories
* using the results of general theory," Astronomy and Astrophysics
* 157, 59070 (1986).
*
* See also:
* P. Bretagnon and G. Francou, "Planetary theories in rectangular
* and spherical variables. VSOP87 solutions," Astronomy and
* Astrophysics 202, 309-315 (1988).
*
* Laskar's expansions are said by Bretagnon and Francou
* to have "a precision of about 1" over 10000 years before
* and after J2000.0 in so far as the precession constants p^0_A
* and epsilon^0_A are perfectly known."
*
* Bretagnon and Francou's expansions for the node and inclination
* of the ecliptic were derived from Laskar's data but were truncated
* after the term in T**6. I have recomputed these expansions from
* Laskar's data, retaining powers up to T**10 in the result.
*
* The following table indicates the differences between the result
* of the IAU formula and Laskar's formula using four different test
* vectors, checking at J2000 plus and minus the indicated number
* of years.
*
* Years Arc
* from J2000 Seconds
* ---------- -------
* 0 0
* 100 .006
* 200 .006
* 500 .015
* 1000 .28
* 2000 6.4
* 3000 38.
* 10000 9400.
*/
/***** modification history
*
* $Log: precess.c,v $
* Revision 1.3 1993/04/21 21:43:02 craig
* Changed the path of the satellite.h include.
* Changed ecnsts to pcnsts.
*
* Revision 1.2 1993/04/21 15:32:40 craig
* First working version. Ran through indent and coverted to ansi.
* Added hooks for working with the satellite programs.
*
*
*/
/***** include files *****/
#include <math.h>
#include "aaproto.h"
#include "satellite.h"
/***** global variables *****/
extern struct MCONSTANTS mcnsts;
extern struct PCONSTANTS pcnsts;
/* from epsiln.c */
extern double coseps, sineps;
/***** local global variables *****/
/* Precession coefficients taken from Laskar's paper: */
static double pAcof[] = {
-8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
-0.235316, 0.07732, 111.1971, 50290.966};
/* Node and inclination of the earth's orbit computed from
* Laskar's data as done in Bretagnon and Francou's paper:
*/
static double nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10,
-3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
-0.042078604317, 3.052112654975};
static double inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
-5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
0.002278495537, 0.0};
/****************************************************************/
/* precess */
/* */
/* Subroutine arguments: */
/* */
/* R = rectangular equatorial coordinate vector to be precessed.*/
/* The result is written back into the input vector. */
/* J = Julian date */
/* direction = */
/* Precess from J to J2000: direction = 1 */
/* Precess from J2000 to J: direction = -1 */
/* Note that if you want to precess from J1 to J2, you would */
/* first go from J1 to J2000, then call the program again */
/* to go from J2000 to J2. */
/****************************************************************/
int precess (double *R, double J, int direction)
{
int i;
double sinth, costh, sinZ, cosZ, sinz, cosz;
double A, B, T, Z, z, TH, pA, W;
double x[3];
double *p;
if (J == pcnsts.J2000)
{
return (0);
}
/* Each precession angle is specified by a polynomial in T = Julian
centuries from J2000.0. See AA page B18. */
T = (J - pcnsts.J2000) / pcnsts.dapcen;
if (fabs (T) > 2.0)
{
goto laskar;
}
Z = ((0.017998 * T + 0.30188) * T + 2306.2181) * T * mcnsts.sec2ra;
z = ((0.018203 * T + 1.09468) * T + 2306.2181) * T * mcnsts.sec2ra;
TH = ((-0.041833 * T - 0.42665) * T + 2004.3109) * T * mcnsts.sec2ra;
sinth = sin (TH);
costh = cos (TH);
sinZ = sin (Z);
cosZ = cos (Z);
sinz = sin (z);
cosz = cos (z);
A = cosZ * costh;
B = sinZ * costh;
if (direction < 0)
{
/* From J2000.0 to J */
x[0] = (A * cosz - sinZ * sinz) * R[0]
- (B * cosz + cosZ * sinz) * R[1]
- sinth * cosz * R[2];
x[1] = (A * sinz + sinZ * cosz) * R[0]
- (B * sinz - cosZ * cosz) * R[1]
- sinth * sinz * R[2];
x[2] = cosZ * sinth * R[0]
- sinZ * sinth * R[1]
+ costh * R[2];
}
else
{
/* From J to J2000.0 */
x[0] = (A * cosz - sinZ * sinz) * R[0]
+ (A * sinz + sinZ * cosz) * R[1]
+ cosZ * sinth * R[2];
x[1] = -(B * cosz + cosZ * sinz) * R[0]
- (B * sinz - cosZ * cosz) * R[1]
- sinZ * sinth * R[2];
x[2] = -sinth * cosz * R[0]
- sinth * sinz * R[1]
+ costh * R[2];
}
goto done;
laskar:
/* Implementation by elementary rotations using Laskar's expansions.
First rotate about the x axis from the initial equator to the
ecliptic. (The input is equatorial.) */
if (direction == 1)
{
epsiln (J); /* To J2000 */
}
else
{
epsiln (pcnsts.J2000); /* From J2000 */
}
x[0] = R[0];
z = coseps * R[1] + sineps * R[2];
x[2] = -sineps * R[1] + coseps * R[2];
x[1] = z;
/* Precession in longitude */
T /= 10.0; /* thousands of years */
p = pAcof;
pA = *p++;
for (i = 0; i < 9; i++)
{
pA = pA * T + *p++;
}
pA *= mcnsts.sec2ra * T;
/* Node of the moving ecliptic on the J2000 ecliptic. */
p = nodecof;
W = *p++;
for (i = 0; i < 10; i++)
{
W = W * T + *p++;
}
/* Rotate about z axis to the node. */
if (direction == 1)
{
z = W + pA;
}
else
{
z = W;
}
B = cos (z);
A = sin (z);
z = B * x[0] + A * x[1];
x[1] = -A * x[0] + B * x[1];
x[0] = z;
/* Rotate about new x axis by the inclination of the moving ecliptic
on the J2000 ecliptic. */
p = inclcof;
z = *p++;
for (i = 0; i < 10; i++)
{
z = z * T + *p++;
}
if (direction == 1)
{
z = -z;
}
B = cos (z);
A = sin (z);
z = B * x[1] + A * x[2];
x[2] = -A * x[1] + B * x[2];
x[1] = z;
/* Rotate about new z axis back from the node. */
if (direction == 1)
{
z = -W;
}
else
{
z = -W - pA;
}
B = cos (z);
A = sin (z);
z = B * x[0] + A * x[1];
x[1] = -A * x[0] + B * x[1];
x[0] = z;
/* Rotate about x axis to final equator. */
if (direction == 1)
{
epsiln (pcnsts.J2000);
}
else
{
epsiln (J);
}
z = coseps * x[1] - sineps * x[2];
x[2] = sineps * x[1] + coseps * x[2];
x[1] = z;
done:
for (i = 0; i < 3; i++)
{
R[i] = x[i];
}
return (0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.