/* Expansions for the geocentric ecliptic longitude, latitude,
 * and equatorial horizontal parallax of the Moon referred to
 * the mean equinox of date.
 *
 * Coefficients are from Jean Meeus, "Astronomical Formulae
 * for Calculators," 3rd ed., Willman-Bell, 1986.
 * Accuracy estimated by Meeus is 10" longitude, 3" latitude,
 * and 0.2" parallax.  A check against the 1986 Almanac at
 * 35 different dates revealed peak error of 0.5s R.A.
 * and 4" Declination, in agreement.
 *
 * Corrections for nutation and light time are included.
 * Note that the difference between Ephemeris and Universal
 * time is significant here, since the Moon moves in right
 * ascension at the rate of about 2s/minute.
 *
 * - S. L. Moshier, October, 1987
 */

#include "kep.h"

/* Rate of motion in R.A. and Dec., computed by this program.
 */
extern double dradt, ddecdt;

/* Obliquity of the ecliptic, computed by epsiln():
 */
extern double coseps, sineps;

/* Answers posted by angles():
 */
extern double pq, qe, ep;

/* Nutation, computed by nutlo():
 */
extern double nutl, nuto;

/* The following times are set up by update() and refer
 * to the same instant.  The distinction between them
 * is required by altaz().
 */
extern double TDT; /* Ephemeris time */
extern double UT;  /* Universal time */

#define DEBUG 0
static double l = 0.0;		/* Moon's ecliptic longitude */
static double B = 0.0;		/* Ecliptic latitude */
static double p = 0.0;		/* Parallax */
static double ra = 0.0;		/* Right Ascension */
static double dec = 0.0;	/* Declination */

/* Beware of atan2()!
 */
double floor(), sin(), cos(), asin();

/* equatorial radius of the earth, in au
 */
#define Rearth 4.263523e-5


/* Program begins:
 */
domoon()
{
int i;
double ra0, dec0, r;
double x, y, z, lon0;
double pp[3], qq[3];
double acos();

/* Compute obliquity of the ecliptic, coseps, and sineps
 * (see kep.c).
 */
epsiln(TDT);

/* Run the orbit calculation twice, at two different times,
 * in order to find the rate of change of R.A. and Dec.
 */

/* Calculate for 0.001 day ago
 */
i = prtflg;
prtflg = 0; /* disable display */
moonll(TDT-0.001);
lonlat( rearth, TDT, eapolar ); /* precess earth to date */
prtflg = i; /* reenable display */
ra0 = ra;
dec0 = dec;
lon0 = l;

/* Calculate for present instant.
 */
moonll(TDT);

/* The rates of change.  These are used by altaz() to
 * correct the time of rising, transit, and setting.
 */
dradt = 1000.0*(ra-ra0);
ddecdt = 1000.0*(dec-dec0);

/* Post the ecliptic longitude and latitude, in radians,
 * and the radius in au.
 */
obpolar[0] = l;
obpolar[1] = B;
r = 1.0/sin(p); /* distance in earth-radii */
ra0 = 4.263523e-5*r; /* factor is radius of earth in au */
obpolar[2] = ra0;

/* Rate of change in longitude, degrees per day
 * used for phase of the moon
 */
lon0 = 1000.0*RTD*(l - lon0);

/* convert to ecliptic rectangular coordinates */

z = ra0 * cos(B);
x = z * cos(l);
y = z * sin(l);
z = ra0 * sin(B);
/* convert to equatorial coordinates */
pp[0] = x;
pp[1] = y * coseps - z * sineps;
pp[2] = y * sineps + z * coseps;

/* Find sun-moon-earth angles */
precess( rearth, TDT, -1 );
for( i=0; i<3; i++ )
	qq[i] = rearth[i] + pp[i];
angles( pp, qq, rearth );

/* Display answers
 */
if( prtflg )
	{
	printf( "Apparent geocentric longitude %.3lf deg", RTD*(l+nutl) );
	printf( "   latitude %.3lf deg\n", RTD*B );
	printf( "Distance %.5lf Earth-radii\n", r );
	printf( "Horizontal parallax" );
	dms( p );
	printf( "Semidiameter" );
	x = 0.272453 * p  +  0.0799/RTS; /* AA page L6 */
	dms( x );

	x = RTD * acos(-ep);
	printf( "\nElongation from sun %.2lf deg,", x );
	x = 0.5 * (1.0 + pq);
	printf( "  Illuminated fraction %.2lf\n", x );

/* Find phase of the Moon by comparing Moon's longitude
 * with Earth's longitude.
 *
 * The number of days before or past indicated phase is
 * estimated by assuming the true longitudes change linearly
 * with time.  These rates are estimated for the date, but
 * do not stay constant.  The error can exceed 0.15 day in 4 days.
 */
/* Apparent longitude of sun.
 * Moon's longitude was already corrected for light time.
 */
	y = eapolar[0] - 20.496/(RTS*eapolar[2]);
	x = obpolar[0] - y;
	x = modtp( x ) * RTD;	/* difference in longitude */
	i = x/90;		/* number of quarters */
	x = (x - i*90.0);	/* phase angle mod 90 degrees */

/* days per degree of phase angle */
	z = 1.0/(lon0 - (0.9856/eapolar[2]));

	if( x > 45.0 )
		{
		y = -(x - 90.0)*z;
		printf( "Phase %.1lf days before ", y );
		i = (i+1) & 3;
		}
	else
		{
		y = x*z;
		printf( "Phase %.1lf days past ", y );
		}


	switch(i)
		{
		case 0: printf( "Full Moon\n" ); break;
		case 1: printf( "Third Quarter\n" ); break;
		case 2: printf( "New Moon\n" ); break;
		case 3: printf( "First Quarter\n" ); break;
		}
	} /* if prtflg */

printf( "    Apparent:  R.A." );
hms(ra);
printf( "Declination" );
dms(dec);
printf( "\n" );

/* Compute and display topocentric position (altaz.c)
 */
pp[0] = ra;
pp[1] = dec;
pp[2] = r * Rearth;
altaz( pp, UT );
}



/* Orbit calculation begins.
 */
static double e = 0.0;
static double ee = 0.0;
static double f = 0.0;
static double LP = 0.0;
static double M = 0.0;
static double MP = 0.0;
static double D = 0.0;
static double NF = 0.0;
static double OM = 0.0;
static double T = 0.0;
static double om2 = 0.0;
static double smp = 0.0;
static double cmp = 0.0;
static double s2mp = 0.0;
static double c2mp = 0.0;
static double sd = 0.0;
static double cd = 0.0;
static double s2d = 0.0;
static double c2d = 0.0;
static double sf = 0.0;
static double cf = 0.0;
static double s2f = 0.0;
static double c2f = 0.0;
static double sm = 0.0;
static double cm = 0.0;
static double s2m = 0.0;
static double c2m = 0.0;
static double sx = 0.0;
static double cx = 0.0;
static double sy = 0.0;
static double cy = 0.0;
static double sz = 0.0;
static double cz = 0.0;


/* Calculate latitude, longitude, and horizontal parallax
 * of the Moon at Julian date J (time scale is TDT).
 */
moonll(J)
double J;
{

T = (J-J1900)/36525.0;
/* Program is broken up to get it thru micro compiler */
moon1();
moon2();
moon3();
moon4();
}



moon1()
{
double a, b;

/* Mean longitude of moon */
LP = ((1.9e-6*T - 0.001133)*T + 481267.8831)*T + 270.434164;
/* Mean anomaly of sun */
M = (( -3.3e-6*T - 1.50e-4)*T + 35999.0498)*T + 358.475833;
/* Mean anomaly of moon */
MP = ((1.44e-5*T + 0.009192)*T + 477198.8491)*T + 296.104608;
/* Mean elongation of moon */
D = ((1.9e-6*T - 0.001436)*T + 445267.1142)*T + 350.737486;
/* Mean distance of moon from its ascending node */
NF = (( -3.e-7*T - 0.003211)*T + 483202.0251)*T + 11.250889;
/* Longitude of ascending node of moon */
OM = (( 2.2e-6*T + 0.002078)*T - 1934.1420)*T + 259.183275;

f = 51.2 + 20.2*T;
f = sin(DTR*f);
LP += 0.000233*f;
M  -= 0.001778*f;
MP += 0.000817*f;
D  += 0.002011*f;

f = ( -0.0091731*T + 132.870)*T + 346.560;
f = 0.003964*sin(DTR*f);
LP += f;
MP += f;
D  += f;
NF  += f;

f = sin(DTR*OM);
LP += 0.001964*f;
MP += 0.002541*f;
D  += 0.001964*f;
NF  -= 0.024691*f; 

f = (-2.30*T + 275.05 + OM)*DTR;
NF -= 0.004328*sin(f);
om2 = 0.0000754*cos(f);

e = ( -7.52e-6*T - 0.002495)*T + 1.0;
ee = e * e;

l = LP;

/* convert all to radians */
M  *= DTR;
MP *= DTR;
D  *= DTR;
NF  *= DTR;

sm = sin(M);
cm = cos(M);

/* We will save a lot of computer time by using
 * trigonometric identities instead of
 * the sin() and cos() functions.
 */
s2m = 2.0*sm*cm;		/* sin(2M) */
c2m = cm*cm - sm*sm;		/* cos(2M) */
smp = sin(MP);
cmp = cos(MP);
s2mp = 2.0*smp*cmp;		/* sin(2MP) */
c2mp = cmp*cmp - smp*smp;	/* cos(2MP) */
sd = sin(D);
cd = cos(D);
s2d = 2.0*sd*cd;	/* sin(2D) */
c2d = cd*cd - sd*sd;	/* cos(2D) */
sf = sin(NF);
cf = cos(NF);
s2f = 2.0*sf*cf;	/* sin(2F) */
c2f = cf*cf - sf*sf;	/* cos(2F) */
sx = s2d*cmp - c2d*smp;	/* sin(2D - MP) */
cx = c2d*cmp + s2d*smp;	/* cos(2D - MP) */
sy = sx*cm - cx*sm;	/* sin(2D-MP-M) */
cy = cx*cm + sx*sm;	/* cos(2D-MP-M) */
sz = s2d*c2mp - c2d*s2mp;	/* sin(2D-2MP) */
cz = c2d*c2mp + s2d*s2mp;	/* cos(2D-2MP) */
l = l +  6.288750*smp
	+1.274018*sx
	+0.658309*s2d
	+0.213616*s2mp
	-0.185596*sm*e;
l = l -  0.114336*s2f
	+0.058793*sz
	+0.057212*sy*e
	-0.034718*sd
	-0.002079*s2m*ee;

a = smp*cf;
b = cmp*sf;
B =	 5.128189*sf
	+0.280606*(a+b)		/* sin(MP+F) */
	+0.277693*(a-b);	/* sin(MP-F) */

p =	 0.950724
	+0.051818*cmp
	+0.009531*cx
	+0.007843*c2d
	+0.002824*c2mp;
p = p +  0.000401*cy*e
	-0.000271*cd
	-0.000084*cz
	-0.000111*cm*e;
a = sz*cm;
b = cz*sm;
l += 0.000693*e*(a+b);		/* sin(2D-2MP+M) */
l += 0.002396*e*(a-b);		/* sin(2D-2MP-M) */
a = sz*cf;
b = cz*sf;
B -= 0.000450*(a+b);		/* sin(2D-2MP+F) */
B += 0.004323*(a-b);		/* sin(2D-2MP-F) */
}



moon2()
{
double a, b;

sz = s2d*cf - c2d*sf;	/* sin(2D-F) */
B +=	0.173238*sz;
cz = c2d*cf + s2d*sf;	/* cos(2D-F) */
f = sz*cm + cz*sm;	/* sin(2D-F+M) */
B -= 0.003372*f*e;
a = sx*cf;
b = cx*sf;
B +=	0.055413*(a+b);	/* sin(2D-MP+F) */
B +=	0.046272*(a-b);	/* sin(2D-MP-F) */
B += 0.032573*(s2d*cf + c2d*sf); /* sin(2D+F) */
B += 0.017198*(s2mp*cf + c2mp*sf);	/* sin(2MP+F) */
a = sy*cf;
b = cy*sf;
B +=	0.002472*e*(a+b);	/* sin(2D-MP-M+F) */
B +=	0.002072*e*(a-b);	/* sin(2D-MP-M-F) */

sz = sx*cm + cx*sm;	/* sin(2D-MP+M) */
l -= 0.007910*sz*e;
cz = cx*cm - sx*sm;	/* cos(2D-MP+M) */
p -= 0.000063*cz*e;
B -= 0.000367*e*(sz*cf + cz*sf); /* sin(2D-MP+M+F) */

l -= 0.002602*(sx*c2f + cx*s2f); /* sin(2D-MP+2F) */
a = sx*c2m;
b = cx*s2m;
l -= 0.000704*ee*(a+b);	/* sin(2D-MP+2M) */
l += 0.002059*ee*(a-b);	/* sin(2D-MP-2M) */
p -= 0.000023*(cx*c2f + sx*s2f); /* cos(2D-MP-2F) */

sx = s2d*cmp + c2d*smp;	/* sin(2D+MP) */
l += 0.053320*sx;
cx = c2d*cmp - s2d*smp;	/* cos(2D+MP) */
p += 0.000857*cx;
a = sx*cm;
b = cx*sm;
l -= 0.000811*e*(a+b);	/* sin(2D+MP+M) */
l += 0.004049*e*(a-b);	/* sin(2D+MP-M) */
p += 0.000064*e*(cx*cm + sx*sm); /* cos(2D+MP-M) */

sz = sx*cf - cx*sf;	/* sin(2D+MP-F) */
B += 0.009267*sz;
cz = cx*cf + sx*sf;	/* cos(2D+MP-F) */
B += 0.000492*e*(sz*cm - cz*sm); /* sin(2D+MP-F-M) */

sz = sx*cf + cx*sf;	/* sin(2D+MP+F) */
B += 0.004200*sz;
cz = cx*cf - sx*sf;	/* cos(2D+MP+F) */
B += 0.000317*e*(sz*cm - cz*sm); /* sin(2D+MP+F-M) */
l -= 0.001773*(sx*c2f - cx*s2f); /* sin(2D+MP-2F) */
B += 0.008823*(s2mp*cf - c2mp*sf); /* sin(2MP-F) */

sz = s2d*cm - c2d*sm;	/* sin(2D-M) */
l += 0.045874*sz*e;
cz = c2d*cm + s2d*sm;	/* cos(2D-M) */
p += 0.000533*cz*e;
a = sz*cf;
b = cz*sf;
B += 0.002222*e*(a+b);	/* sin(2D-M+F) */
B += 0.008247*e*(a-b);	/* sin(2D-M-F) */

sy = smp*cm - cmp*sm;	/* sin(MP-M) */
l += 0.041024*sy*e;
cy = cmp*cm + smp*sm;	/* cos(MP-M) */
p += 0.000320*cy*e;
a = sy*cf;
b = cy*sf;
B += 0.001877*e*(a+b);	/* sin(MP-M+F) */
B += 0.001570*e*(a-b);	/* sin(MP-M-F) */

sz = smp*cm + cmp*sm;	/* sin(MP+M) */
l -= 0.030465*sz*e;
cz = cmp*cm - smp*sm;	/* cos(MP+M) */
p -= 0.000264*cz*e;
a = sz*cf;
b = cz*sf;
B -= 0.001481*e*(a+b);	/* sin(MP+M+F) */
B -= 0.001417*e*(a-b);	/* sin(MP+M-F) */

sy = s2d*c2f - c2d*s2f;	/* sin(2D-2F) */
l += 0.015326*sy;
cy = c2d*c2f + s2d*s2f;	/* cos(2D-2F) */
p -= 0.000029*cy;
l += 0.000598*e*(sy*cm - cy*sm); /* sin(2D-2F-M) */

a = s2f*cmp;
b = c2f*smp;
l -= 0.012528*(a+b);	/* sin(2F+MP) */
l -= 0.010980*(a-b);	/* sin(2F-MP) */
p -= 0.000198*(c2f*cmp + s2f*smp); /* cos(2F-MP) */

sy = 2.0*s2d*c2d;	/* sin(4D) */
l += 0.003862*sy;
cy = c2d*c2d - s2d*s2d;	/* cos(4D) */
p += 0.000072*cy;
l += 0.000550*(sy*cmp + cy*smp); /* sin(4D+MP) */

sz = sy*cmp - cy*smp;	/* sin(4D-MP) */
l += 0.010674*sz;
cz = cy*cmp + sy*smp;	/* cos(4D-MP) */
p += 0.000167*cz;
a = sz*cf;
b = cz*sf;
B += 0.000833*(a+b);	/* sin(4D-MP+F) */
B += 0.001828*(a-b);	/* sin(4D-MP-F) */
}

moon3()
{
double a, b;

sz = sy*cm - cy*sm;	/* sin(4D-M) */
l += 0.000521*sz*e;
cz = cy*cm + sy*sm;	/* cos(4D-M) */
l += 0.000761*e*(sz*c2mp - cz*s2mp); /* sin(4D-M-2MP) */
l += 0.001220*e*(sz*cmp - cz*smp); /* sin(4D-M-MP) */
p += 0.000019*e*(cz*cmp + sz*smp); /* cos(4D-M-MP) */

sz = sy*c2mp - cy*s2mp;	/* sin(4D-2MP) */
l += 0.008548*sz;
cz = cy*c2mp + sy*s2mp;	/* cos(4D-2MP) */
p += 0.000103*cz;
B += 0.000670*(sz*cf + cz*sf); /* sin(4D-2MP+F) */
a = sy*cf;
b = cy*sf;
B += 0.001020*(a-b);	/* sin(4D-F) */
B += 0.000331*(a+b);	/* sin(4D+F) */

sy = s2mp*cmp + c2mp*smp;	/* sin(3MP) */
l += 0.010034*sy;
cy = c2mp*cmp - s2mp*smp;	/* cos(3MP) */
p += 0.000173*cy;
l -= 0.003665*(sy*c2d - cy*s2d); /* sin(3MP-2D) */
p -= 0.000033*(cy*c2d + sy*s2d); /* cos(3MP-2D) */
a = sy*cf;
b = cy*sf;
B += 0.001106*(a+b);		/* sin(3MP+F) */
B += 0.000439*(a-b);		/* sin(3MP-F) */

a = 2.0*D - NF - 3.0*MP;
B += 0.000422*sin(a);		/* sin(2D-F-3MP) */

sy = s2f*cf + c2f*sf;		/* sin(3F) */
B -= 0.001750*sy;
cy = c2f*cf - s2f*sf;		/* cos(3F) */

B += 0.000606*(s2d*cy - c2d*sy); /* sin(2D-3F) */
a = smp*cy;
b = cmp*sy;
B += 0.000781*(a-b);		/* sin(MP-3F) */
B -= 0.000283*(a+b);		/* sin(MP+3F) */
a = sf*cm;
b = cf*sm;
B -= 0.001803*e*(a+b);		/* sin(F+M) */
B += 0.001350*e*(a-b);		/* sin(F-M) */
a = sf*cd;
b = cf*sd;
B -= 0.001487*(a+b);		/* sin(F+D) */
B += 0.001330*(a-b);		/* sin(F-D) */
a = smp*cd;
b = cmp*sd;
l += 0.005162*(a-b);		/* sin(MP-D) */
l -= 0.002349*(a+b);		/* sin(MP+D) */
p -= 0.000030*(cmp*cd - smp*sd); /* cos(MP+D) */
l += 0.005000*e*(sm*cd + cm*sd); /* sin(M+D) */
p += 0.000041*e*(cm*cd - sm*sd); /* cos(M+D) */
a = s2mp*cm;
b = c2mp*sm;
l -= 0.002125*e*(a+b);		/* sin(2MP+M) */
l += 0.002695*e*(a-b);		/* sin(2MP-M) */
a = c2mp*cm;
b = s2mp*sm;
p -= 0.000029*e*(a-b);		/* cos(2MP+M) */
p += 0.000035*e*(a+b);		/* cos(2MP-M) */

sz = s2d*c2m - c2d*s2m;		/* sin(2D-2M) */
l += 0.002249*sz*ee;
cz = c2d*c2m + s2d*s2m;		/* cos(2D-2M) */
p += 0.000026*cz*ee;
B += 0.000306*ee*(sz*cf - cz*sf); /* sin(2D-2M-F) */
l -= 0.001595*(s2d*c2f + c2d*s2f); /* sin(2D+2F) */
l -= 0.001110*(s2mp*c2f + c2mp*s2f); /* sin(2MP+2F) */

sz = s2d*c2mp + c2d*s2mp;	/* sin(2D+2MP) */
l += 0.003996*sz;
cz = c2d*c2mp - s2d*s2mp;	/* cos(2D+2MP) */
p += 0.000079*cz;
a = sz*cf;
b = cz*sf;
B += 0.000423*(a+b);		/* sin(2D+2MP+F) */
B += 0.000597*(a-b);		/* sin(2D+2MP-F) */

sy = s2d*cm + c2d*sm;		/* sin(2D+M) */
l -= 0.006783*sy*e;
cy = c2d*cm - s2d*sm;		/* cos(2D+M) */
p -= 0.000083*cy*e;
B -= 0.000353*e*(sy*cf + cy*sf); /* sin(2D+M+F) */
l += 0.000717*ee*(smp*c2m - cmp*s2m); /* sin(MP-2M) */
l += 0.001076*s2mp*c2mp; /* sin(4MP) */
l += 0.000486*(s2mp*cd - c2mp*sd); /* sin(2MP-D) */

l += 0.000892*sin(MP-3.0*D);

l -= 360.0 * floor( l/360.0 );

B = B * (1.0 - 0.0004664*cos(OM*DTR) - om2 );

l *= DTR;
B *= DTR;
p *= DTR;
}


/* Convert ecliptic polar coordinates
 * to R.A. and Dec.
 */
moon4()
{
double cosB, sinB, cosL, sinL;
double pp[3];


/* Light time correction to longitude,
 * about 0.7".
 */
l -= DTR * 0.0118 * sin(p);

/* Subtract 0.8 arcsec per century squared from the longitude
 * to modernize the estimated tidal acceleration of the Moon's orbit.
 */
l -= 3.88e-6 * T * T;

cosB = cos(B);
sinB = sin(B);
cosL = cos(l);
sinL = sin(l);

pp[0] = cosB*cosL;
pp[1] = coseps*cosB*sinL - sineps*sinB;
pp[2] = sineps*cosB*sinL + coseps*sinB;

/* Correct for nutation at date TDT.
 */
nutate( TDT, pp );


ra = zatan2(pp[0],pp[1]);
dec = asin(pp[2]);
}