/* MoonRise.c
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA 02111-1307, USA.
 *
 * Original (C) Mike Henderson <mghenderson@lanl.gov>.
 *
 * I've converted Mike Henderson's original file to glib types, 
 * removed unused variables, and piped the whole thing through indent.
 *
 * (C) 1999-2003 josh buhl <jbuhl@users.sourceforge.net> 
 *
 */

#include "Moon.h"
#include "MoonRise.h"

static gint
Interp(gdouble ym, gdouble y0, gdouble yp, gdouble * xe, gdouble * ye,
       gdouble * z1, gdouble * z2, gint * nz)
{

	gdouble a, b, c, d, dx;

	*nz = 0;
	a = 0.5 * (ym + yp) - y0;
	b = 0.5 * (yp - ym);
	c = y0;
	*xe = -b / (2.0 * a);
	*ye = (a * (*xe) + b) * (*xe) + c;
	d = b * b - 4.0 * a * c;

	if (d >= 0) {
		dx = 0.5 * sqrt(d) / fabs(a);
		*z1 = *xe - dx;
		*z2 = *xe + dx;
		if (fabs(*z1) <= 1.0)
			*nz += 1;
		if (fabs(*z2) <= 1.0)
			*nz += 1;
		if (*z1 < -1.0)
			*z1 = *z2;
	}

	return (0);
}


static gdouble 
SinH(CTrans * c, gdouble UT)
{

	gdouble TU, TU2, TU3, LambdaMoon, BetaMoon, R, AGE, frac(), jd();
	gdouble RA_Moon, DEC_Moon, gmst, lmst, Tau, epsilon, Moon();
	gdouble angle2pi();

	TU = (jd(c->year, c->month, c->day, UT) - 2451545.0) / 36525.0;

	/* this is more accurate, but expensive. */
	TU2 = TU * TU;
	TU3 = TU2 * TU;
	Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
	LambdaMoon *= RadPerDeg;
	BetaMoon *= RadPerDeg;
	epsilon =
	    (23.43929167 - 0.013004166 * TU - 1.6666667e-7 * TU2 -
	     5.0277777778e-7 * TU3) * RadPerDeg;
	RA_Moon =
	    angle2pi(atan2
		     (sin(LambdaMoon) * cos(epsilon) -
		      tan(BetaMoon) * sin(epsilon), cos(LambdaMoon)));
	DEC_Moon =
	    asin(sin(BetaMoon) * cos(epsilon) +
		 cos(BetaMoon) * sin(epsilon) * sin(LambdaMoon));


	/* This is less accurate, but computationally less expensive */
/*      MiniMoon(TU, &RA_Moon, &DEC_Moon); */
/*      RA_Moon *= 15.0*RadPerDeg; */
/*      DEC_Moon *= RadPerDeg; */

	/*
	 *  Compute Greenwich Mean Sidereal Time (gmst)
	 */
	UT = 24.0 * frac(UT / 24.0);
	/* this is for the ephemeris meridian??? 
	   gmst = 6.697374558 + 1.0027379093*UT + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
	 */
	gmst =
	    UT + 6.697374558 + (8640184.812866 +
				(0.093104 -
				 6.2e-6 * TU) * TU) * TU / 3600.0;
	lmst = 24.0 * frac((gmst - c->Glon / 15.0) / 24.0);

	Tau = 15.0 * lmst * RadPerDeg - RA_Moon;

	return (c->SinGlat * sin(DEC_Moon) +
		c->CosGlat * cos(DEC_Moon) * cos(Tau));
}


void MoonRise(CTrans * c, gdouble * UTRise, gdouble * UTSet)
{

	gdouble UT, ym, y0, yp, SinH0;
	gdouble xe, ye, z1, z2;
	gint Rise, Set, nz, TimeZone;

	SinH0 = sin(8.0 / 60.0 * RadPerDeg);

	/* report rise and set times in GMT */
	/*      TimeZone = c->UT - c->LST; */
	TimeZone = 0;

	UT = 1.0 + TimeZone;
	UT = 1.0;
	*UTRise = -999.0;
	*UTSet = -999.0;
	Rise = Set = 0;

	ym = SinH(c, UT - 1.0) - SinH0;

	while ((UT <= 24.0 + TimeZone)) {

		y0 = SinH(c, UT) - SinH0;
		yp = SinH(c, UT + 1.0) - SinH0;

		Interp(ym, y0, yp, &xe, &ye, &z1, &z2, &nz);

		switch (nz) {

		case 0:
			break;
		case 1:
			if (ym < 0.0) {
				*UTRise = UT + z1;
				Rise = 1;
			} else {
				*UTSet = UT + z1;
				Set = 1;
			}
			break;
		case 2:
			if (ye < 0.0) {
				*UTRise = UT + z2;
				*UTSet = UT + z1;
			} else {
				*UTRise = UT + z1;
				*UTSet = UT + z2;
			}
			Rise = 1;
			Set = 1;
			break;
		}
		ym = yp;
		UT += 2.0;

	}

	if (Rise) {
		*UTRise -= TimeZone;
		*UTRise = hour24(*UTRise);
	} else {
		*UTRise = -999.0;
	}

	if (Set) {
		*UTSet -= TimeZone;
		*UTSet = hour24(*UTSet);
	} else {
		*UTSet = -999.0;
	}

}