/* Precession of the equinox and ecliptic
 * from epoch Julian date J to or from J2000.0
 *
 * Program by Steve Moshier.
 *
 *
 * 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.
 */


/* 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 };

extern double J2000; /* = 2451545.0, 2000 January 1.5 */
extern double STR; /* = 4.8481368110953599359e-6 radians per arc second */
extern double coseps, sineps; /* see epsiln.c */

/* 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.
 */

precess( R, J, direction )
double R[], J;
int direction;
{
double sinth, costh, sinZ, cosZ, sinz, cosz;
double A, B, T, Z, z, TH, pA, W;
double x[3];
double *p;
int i;
double sin(), cos(), fabs();

if( J == J2000 )
	return;
/* Each precession angle is specified by a polynomial in
 * T = Julian centuries from J2000.0.  See AA page B18.
 */
T = (J - J2000)/36525.0;

if( fabs(T) > 2.0 )
	goto laskar;

Z =  (( 0.017998*T + 0.30188)*T + 2306.2181)*T*STR;
z =  (( 0.018203*T + 1.09468)*T + 2306.2181)*T*STR;
TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*STR;

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( 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 *= STR * 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( 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];
}