/* One-term nutation program uses only the 18.6 year term.
 *
 *
 * References:
 * "Summary of 1980 IAU Theory of Nutation (Final Report of the
 * IAU Working Group on Nutation)", P. K. Seidelmann et al., in
 * Transactions of the IAU Vol. XVIII A, Reports on Astronomy,
 * P. A. Wayman, ed.; D. Reidel Pub. Co., 1982.
 *
 * "Nutation and the Earth's Rotation",
 * I.A.U. Symposium No. 78, May, 1977, page 256.
 * I.A.U., 1980.
 *
 * Woolard, E.W., "A redevelopment of the theory of nutation",
 * The Astronomical Journal, 58, 1-3 (1953).
 */

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

/* The answers are posted here by nutlo():
 */
DOUBLE jdnut = -1.0;	/* time to which the nutation applies */
DOUBLE nutl = 0.0;	/* nutation in longitude (radians) */
DOUBLE nuto = 0.0;	/* nutation in obliquity (radians) */
extern DOUBLE eps, coseps, sineps, DTR, J1900, J2000;
DOUBLE FLOOR();
#define mod360(x) (x - ((DOUBLE )360.0)*FLOOR(x/((DOUBLE )360.0)))

nutlo(J)
DOUBLE J;
{
DOUBLE f, g, w, T;
DOUBLE OM;
DOUBLE cu, su, cv, sv, cw, sw;
DOUBLE C, D;
int i, j, k, k1, m;
short *p;
DOUBLE SIN(), COS();

if( jdnut == J )
	return(0);
jdnut = J;

#if 1
T = (J-J1900)/36525.0;
OM = ((0.0080*T + 7.48)*T - 6962911.23)*T + 933059.79;
OM = DTR * mod360(OM/3600.0);

/* Woolard series
C = (-0.01737*T - 17.2327)*SIN(OM);
D = ( 0.00091*T +  9.2100)*COS(OM);
 */
/* 1977 IAU series */
C = (-0.01737*T - 17.186)*SIN(OM);
D = ( 0.00091*T +  9.20159)*COS(OM);

#endif

#if 0
/* Julian centuries from 2000 January 1.5,
 * barycentric dynamical time
 */
T = (J-J2000)/36525.0;

/* Fundamental arguments in the FK5 reference system.
/* The coefficients, originally given to 0.001",
 * are converted here to degrees.
 */

/* longitude of the mean ascending node of the lunar orbit
 * on the ecliptic, measured from the mean equinox of date
 */
OM = ((2.22e-6*T + 0.00207833)*T - 1934.1362608)*T + 125.0445222;
OM = DTR * mod360(OM);

C = (-0.01742*T - 17.1996)*SIN(OM);
D = ( 0.00089*T +  9.2025)*COS(OM);
#endif

/*
printf( "nutation: in longitude %.3f\", in obliquity %.3f\"\n", C, D );
*/
/* Save answers, expressed in radians */
nutl = DTR * C / 3600.0;
nuto = DTR * D / 3600.0;
return(1);
}



/* Nutation -- AA page B20
 * using nutation in longitude and obliquity from nutlo()
 * and obliquity of the ecliptic from epsiln()
 * both calculated for Julian date J.
 *
 * p[] = equatorial rectangular position vector of object for
 * mean ecliptic and equinox of date.
 */
static DOUBLE N00, N01, N02, N10, N11, N12, N20, N21, N22;

nutate( J, p )
DOUBLE J;
DOUBLE p[];
{
DOUBLE ce, se, cl, sl, sino, f;
DOUBLE dp[3], p1[3];
int i;
DOUBLE SIN(), COS();

if( J == jdnut )
	goto matmpy;
nutlo(J); /* be sure we calculated nutl and nuto */
epsiln(J); /* and also the obliquity of date */

f = eps + nuto;
ce = COS( f );
se = SIN( f );
sino = SIN(nuto);
cl = COS( nutl );
sl = SIN( nutl );

/* Apply adjustment
 * to equatorial rectangular coordinates of object.
 *
 * This is a composite of three rotations: rotate about x axis
 * to ecliptic of date; rotate about new z axis by the nutation
 * in longitude; rotate about new x axis back to equator of date
 * plus nutation in obliquity.
 */
N00 =   cl;
N01 =  -sl*coseps;
N02 =  -sl*sineps;

N10 = sl*ce;
N11 = cl*coseps*ce + sineps*se;
N12 = -( sino + (1.0-cl)*sineps*se );

N20 = sl*se;
N21 = sino + (cl-1.0)*se*coseps;
N22 = cl*sineps*se + coseps*ce;

matmpy:

p1[0] = N00*p[0] + N01*p[1] + N02*p[2];
p1[1] = N10*p[0] + N11*p[1] + N12*p[2];
p1[2] = N20*p[0] + N21*p[1] + N22*p[2];

/*
for( i=0; i<3; i++ )
	dp[i] = p1[i] - p[i];
showcor( "nutation", p, dp );
*/

for( i=0; i<3; i++ )
	p[i] = p1[i];
}


/* Inverse transformation using previously
 * computed matrix elements.
 */
invnut( p1 )
DOUBLE p1[];
{
DOUBLE p[3];
int i;

p[0] = N00*p1[0] + N10*p1[1] + N20*p1[2];
p[1] = N01*p1[0] + N11*p1[1] + N21*p1[2];
p[2] = N02*p1[0] + N12*p1[1] + N22*p1[2];

for( i=0; i<3; i++ )
	p1[i] = p[i];
}