/* Precession of the equinox and ecliptic
 * from epoch Julian date J to or from J2000.0 (DE200) or B1950 (DE102)
 * Several versions of precession formula are provided.
 */

#include "prec.h"
#include "ssystem.h"

extern DOUBLE J2000; /* = 2451545.0, 2000 January 1.5 */
extern DOUBLE B1950;
extern DOUBLE STR; /* = 4.8481368110953599359e-6 radians per arc second */

static DOUBLE jdprec = -1.1;
static DOUBLE P00, P01, P02, P10, P11, P12, P20, P21, P22;
DOUBLE SIN(), COS(), FABS();

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

void precess( R, J, direction )
DOUBLE R[], J;
int direction;
{
DOUBLE sinth, costh, sinZ, cosZ, sinz, cosz;
DOUBLE A, B, T, Z, z, TH;

#if DE200
if( J == J2000 )
	return;
#else
if( J == B1950 )
	return;
#endif

if( J == jdprec )
	goto matmul;

jdprec = J;

/* Each precession angle is specified by a polynomial in
 * time.  See AA page B18.
 */

#if DE102
T = (J - B1950)/((DOUBLE )36524.2199);
/* Newcomb's formula for B1950 */
Z =  (( 0.018*T + 0.302)*T + 2304.948)*T*STR;
z =  (( 0.019*T + 1.093)*T + 2304.948)*T*STR;
TH = ((-0.042*T - 0.426)*T + 2004.255)*T*STR;
#endif
#if DE118
/* IAU formula for B1950
T = (J - B1950)/((DOUBLE )36525.0);
Z =  (( 0.017998*T + 0.302053)*T + 2305.51978)*T*STR;
z =  (( 0.018203*T + 1.094647)*T + 2305.51978)*T*STR;
TH = ((-0.041833*T - 0.426542)*T + 2004.73750)*T*STR;
 */
/* Newcomb's formula for B1950 */
T = (J - B1950)/((DOUBLE )36524.2199);
Z =  (( 0.018*T + 0.302)*T + 2304.948)*T*STR;
z =  (( 0.019*T + 1.093)*T + 2304.948)*T*STR;
TH = ((-0.042*T - 0.426)*T + 2004.255)*T*STR;
#endif
#if DE200
/* IAU formula for J2000 */
T = (J - J2000)/((DOUBLE )36525.0);
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;
/* Newcomb's formula for J2000
T = (J - J2000)/((DOUBLE )36524.2199);
Z =  (( 0.018*T + 0.302)*T + 2305.646)*T*STR;
z =  (( 0.019*T + 1.093)*T + 2305.646)*T*STR;
TH = ((-0.042*T - 0.426)*T + 2003.829)*T*STR;
 */
#endif

sinth = SIN(TH);
costh = COS(TH);
sinZ = SIN(Z);
cosZ = COS(Z);
sinz = SIN(z);
cosz = COS(z);
A = cosZ*costh;
B = sinZ*costh;

P00 = A*cosz - sinZ*sinz;
P01 = -(B*cosz + cosZ*sinz);
P02 = -sinth*cosz;
P10 = A*sinz + sinZ*cosz;
P11 = -(B*sinz - cosZ*cosz);
P12 = -sinth*sinz;
P20 = cosZ*sinth;
P21 = -sinZ*sinth;
P22 = costh;

matmul:

if( direction < 0 )
	{ /* From basic epoch to J */
	A = P00*R[0] + P01*R[1] + P02*R[2];
	B = P10*R[0] + P11*R[1] + P12*R[2];
	R[2] = P20*R[0] + P21*R[1] + P22*R[2];
	}
else
	{ /* From J to basic epoch */
	A = P00*R[0] + P10*R[1] + P20*R[2];
	B = P01*R[0] + P11*R[1] + P21*R[2];
	R[2] = P02*R[0] + P12*R[1] + P22*R[2];
	}	
R[1] = B;
R[0] = A;
}