/* this file includes routines to calculate the positions of Sun and Moon,
	   from Duffett-Smith, "Practical Astronomy with your Pocket Calculator." */

#include <stdio.h>
#include <math.h>

#ifdef VERBOSE
extern char *progname;
#endif

	/* find the position, RA and Dec in epoch J2000.0, for the Sun at the
	   given Julian Date. Actually, it finds the position in 1990.0 
	   coords, and then precesses those to J2000.0. this is only a 
	   low-precision routine, giving the position accurate to a few 
	   arcminutes */


void sun_radec( /* double jd, *ra, *dec */ );
void moon_radec( /* double jd, *ra, *dec */ );
static void sun_mlam( /* double jd, *m, *lam */ );
static void ecl_to_radec( /* double longitude, latitude, *ra, *dec */ );
int moon_riseset( /* double jd, *jdr, *jds */ );
static int moon_times( /* double jd, *jdr, *jds */ );
void moon_phase( /* double jd, *phase */ );


#ifndef PI
#define PI 3.14159265359
#endif  PI
#define DEGTORAD(x)    ((x)*(PI/180.0))
#define RADTODEG(x)    ((x)*(180.0/PI))
#define DOWNTO360(x)   while ((x) > 360.0) { x -= 360.0; } 
#define UPTO360(x)     while ((x) <   0.0) { x += 360.0; } 
#define TO360(x)       DOWNTO360(x); UPTO360(x)

#define ELONG 279.403303	/* solar mean ecliptic long, 1990.0, degrees */
#define PLONG 282.768422	/* solar long at perigee, 1990.0, degrees */
#define ECCEN   0.016713	/* earth's orbital eccen, 1990.0 */
#define OBLIQ  23.437991	/* mean obliquity of the ecliptic, 1990.0 */
#define BASE   2447891.5	/* JD of Jan 0.0, 1990 */

	/* these variable are set by the moon_radec() routine, and
	   referenced by the moon rise/set routine */
static double mm;		/* corrected lunar anomaly, degrees, epoch 1990.0 */
static double ec;		/* correction for eq. of center, deg., 1990.0 */
static double llong;	/* true lunar longitude */

void
sun_radec(jd, ra, dec)
double jd, *ra, *dec;
{
	double d, n, m, ec, lam, l, e, y, x, a;

	/* calculate sun's mean anomaly and ecliptic longitude */
	sun_mlam(jd, &m, &lam);
	
	/* now convert ecliptic longitude (lam degrees) and ecliptic 
	   latitude (b=0.0 degrees) to RA and Dec in epoch 1990.0 */
	ecl_to_radec(lam, (double) 0.0, &a, &d);

	/* and finally precess the ra and dec to J2000.0 */
	to_j2000(1990.0, a, d, ra, dec);
}

	/* for the given JD, calculate the Sun's mean anomaly, 'm', and
	   ecliptic longitude, 'lam'. */

static void
sun_mlam(jd, m, lam)
double jd, *m, *lam;
{
	double d, n, ec;

	d = jd - BASE;			/* BASE is JD of Jan 0.0, 1990 */
	n = (360.0/365.242191)*d;
	TO360(n);
	*m = n + ELONG - PLONG;
	if (*m < 0.0)
		*m += 360.0;
	ec = (360.0/PI)*ECCEN*sin(DEGTORAD(*m));
	*lam = n + ec + ELONG;
	if (*lam > 360.0)
		*lam -= 360.0;
}

	/* convert the given ecliptic longitude and latitude (in degrees) 
	   into RA and Dec (decimal degrees), for the epoch 1990.0 */

static void
ecl_to_radec(longitude, latitude, ra, dec)
double longitude, latitude, *ra, *dec;
{
	double l, b, d, y, x, a, e;

	l = DEGTORAD(longitude);
	b = DEGTORAD(latitude);
	e = DEGTORAD(OBLIQ);
	d = sin(b)*cos(e) + cos(b)*sin(e)*sin(l);
	d = asin(d);
	y = sin(l)*cos(e) - tan(b)*sin(e);
	x = cos(l);
	a = atan2(y, x);

	*dec = RADTODEG(d);
	*ra = RADTODEG(a);
}
	
	/* find the position of the moon on the given date.  as before, we 
	   calculate the position for epoch 1990.0 and precess the answer to
	   epoch J2000.0. the output RA and Dec are calculate in decimal
	   degrees and passed in the 'ra' and 'dec' parameters.

	   Note that the mm and ec variables are external to this routine -
	   they are static inside this file, so that other routines can 
	   reference them after this routine sets them.  */

#define M_LONG      318.351648	/* moon's mean long. at epoch 1990.0, degrees */
#define M_PERI       36.340410	/* moon's long of perihelion at epoch, deg. */
#define M_NODE      318.510107	/* moon's long of node at epoch, deg. */
#define M_INCL        5.145396	/* inclination of moon's orbit deg. */
#define M_ECCEN       0.054900	/* eccentricity of moon's orbit */
#define M_AXIS   384401.0		/* semi-major axis of orbit, kilometers */
#define M_SIZE        0.5181	/* angular size of moon at dist M_AXIS, deg. */
#define M_PARA        0.9507	/* parallax at dist M_AXIS, deg. */

void
moon_radec(jd, ra, dec)
double jd, *ra, *dec;
{
	double d, sm, slam, n, ev, ae, a3, a4, v, x, y, a, llam, lbeta;
	
	d = jd - BASE;		/* remember BASE is JD of Jan 0, 1990 */
	
	/* now we calculate the sun's mean anomaly and longitude for this date. */
	sun_mlam(jd, &sm, &slam);
	llong = 13.1763966*d + M_LONG;
	TO360(llong);
	mm = llong - 0.1114041*d - M_PERI;
	TO360(mm);
	n = M_NODE - 0.0529539*d;
	TO360(n);
	ev = 2.0*(llong - slam) - mm;
	ev = 1.2739*sin(DEGTORAD(ev));
	ae = 0.1858*sin(DEGTORAD(sm));
	a3 = 0.37*sin(DEGTORAD(sm));
	mm += ev - ae - a3;
	ec = 6.2886*sin(DEGTORAD(mm));
	a4 = 0.214*sin(2.0*DEGTORAD(mm));
	llong += ev + ec - ae + a4;
	v = 0.6583*(sin(2.0*DEGTORAD(llong - slam)));
	llong += v;							/* true lunar longitude */
	n -= 0.16*sin(DEGTORAD(sm));
	y = sin(DEGTORAD(llong - n))*cos(DEGTORAD(M_INCL));
	x = cos(DEGTORAD(llong - n));
	llam = RADTODEG(atan2(y, x));
	llam += n;
	lbeta = RADTODEG(asin(sin(DEGTORAD(llong - n))*sin(DEGTORAD(M_INCL))));

	/* okay, now we have the lunar ecliptic longitude and latitude in 
	   llam and lbeta, so we convert them to RA and Dec (epoch 1990.0) */
	ecl_to_radec(llam, lbeta, &a, &d);

	/* finally, precess the RA and Dec to J2000.0 */
	to_j2000(1990.0, a, d, ra, dec);
}

	/* calculate the times of moonrise and set which are within twelve hours
	   of the given fractional JD.  

	   Note that the returned times are in Julian Date, NOT UT as most
	   other routines use! 

	   returns 0 if all went well, or -1 if the longitude of the observatory
	   could not be read from the configuration file */

int 
moon_riseset(jd, jdr, jds)
double jd, *jdr, *jds;
{
	double dummy;

#ifdef DEBUG
	printf(" into moon_riseset, JD = %.2lf \n", jd);
#endif
	if (moon_times(jd, jdr, jds) < 0)
		return(-1);
	if (jd - *jdr > 0.5) {
#ifdef DEBUG
		printf("  rise time %.2lf too early - use tomorrow's rise \n", *jdr);
#endif
		if (moon_times(jd + 1.0, jdr, &dummy) < 0)
			return(-1);
	}
	else if (*jdr - jd > 0.5) {
#ifdef DEBUG
		printf("  rise time %.2lf too late - use yesterday's rise \n", *jdr);
#endif
		if (moon_times(jd - 1.0, jdr, &dummy) < 0)
			return(-1);
	}
	if (jd - *jds > 0.5) {
#ifdef DEBUG
		printf("  set time  %.2lf too early - use tomorrow's set \n", *jds);
#endif
		if (moon_times(jd + 1.0, &dummy, jds) < 0)
			return(-1);
	}
	else if (*jds - jd > 0.5) {
#ifdef DEBUG
		printf("  set time  %.2lf too late - use yesterday's set \n", *jds);
#endif
		if (moon_times(jd - 1.0, &dummy, jds) < 0)
			return(-1);
	}

	return(0);
}

	/* calculate the times of moonrise and set for the UT date which 
	   corresponds to the given JD.  uses the methods found in Duffett-Smith.
	   Note that the returned times are in Julian Date, NOT UT as most
	   other routines use! 

	   Note also that I have been unable to get this to work correctly, 
	   in the sense that the times which come out of the calculations always
	   seem to be about a half-hour late.  So I add thirty minutes
	   from the times at the end, just to make them come closer to the
	   actual values.  I'm sure that there's a bug in here somewhere,
	   I just don't know where.

	   returns 0 if all went well, or -1 if the longitude of the observatory
	   could not be read from the configuration file */


#define KLUDGE       20.0		/* time in minutes to add to rise/set */
                          		/*   because I can't get it right otherwise */

static int
moon_times(jd, jdr, jds)
double jd, *jdr, *jds;
{
	double longitude, latitude, jd_date, ra1, dec1, ra2, dec2;
	double utr1, uts1, utr2, uts2, gstr1, gsts1, gstr2, gsts2;
	double utr, uts, gsts, gstr;
	double j, x, y, a, rho, theta, pi, t00, psi, dt;
	double rar, decr, ras, decs;

	if (read_longitude(&longitude) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.moon_riseset: can't get longitude from config file\n",
				progname);
#endif
		return(-1);
	}
	if (read_latitude(&latitude) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.moon_riseset: can't get latitude from config file\n",
				progname);
#endif
		return(-1);
	}

	jd_date = floor(jd + 0.5) - 0.5;

	moon_radec(jd_date, &ra1, &dec1);
	obj_riseset(jd_date, ra1, dec1, &utr1, &uts1);
	ut2gst(jd_date, utr1, &gstr1);
	ut2gst(jd_date, uts1, &gsts1);
	moon_radec(jd_date + 0.5, &ra2, &dec2);
	obj_riseset(jd_date, ra2, dec2, &utr2, &uts2);
	ut2gst(jd_date, utr2, &gstr2);
	ut2gst(jd_date, uts2, &gsts2);

	/* t00 is the Greenwich sidereal time at 0 hours on the 
	   observatory longitude, whatever that means. */
	ut2gst(jd_date, 0.0, &t00);
	x = -longitude*1.002738;
	t00 -= x;
	if (t00 < 0.0)
		t00 += 24.0;

	if (gstr1 < t00) {
		gstr1 += 24.0;
		gstr2 += 24.0;
	}
	if (gsts1 < t00) {
		gsts1 += 24.0;
		gsts2 += 24.0;
	}

	/* interpolate to find the approximation to GST of moonrise, set */
	gstr = (12.03*gstr1 - t00*(gstr2 - gstr1))/(12.03 + gstr1 - gstr2);
	gsts = (12.03*gsts1 - t00*(gsts2 - gsts1))/(12.03 + gsts1 - gsts2);
	
	/* now convert to UT and thence JD, so that we can find the distance
	   to the moon when it rises/sets */
	gst2ut(jd_date, gstr, &utr);
	gst2ut(jd_date, gsts, &uts);

	/* okay, calculate the rise time accurately, including corrections
	   for atmospheric refraction (34 arcmin), parallax and size of moon */
	j = jd_date + utr/24.0;
	moon_radec(j, &rar, &decr);
	x = DEGTORAD(mm) + DEGTORAD(ec);
	rho = (1.0 - M_ECCEN*M_ECCEN)/(1.0 + M_ECCEN*cos(x));
	theta = M_SIZE/rho;
	pi = M_PARA/rho;
	x = -pi + (theta/2.0) + 0.5667;		/* vertical displacement, degrees */
	psi = sin(DEGTORAD(latitude))/cos(DEGTORAD(decr));
	psi = acos(psi);
	y = RADTODEG(asin(sin(DEGTORAD(x))/sin(psi)));	/* hor. displ, degrees */
	dt = (240.0*y)/cos(DEGTORAD(decr));		/* rise/set correction, sec */
	gstr -= dt/3600.0;						/* since GST in hours */

	/* ditto for the set time */
	j = jd_date + uts/24.0;
	moon_radec(j, &ras, &decs);
	x = DEGTORAD(mm) + DEGTORAD(ec);
	rho = (1.0 - M_ECCEN*M_ECCEN)/(1.0 + M_ECCEN*cos(x));
	theta = M_SIZE/rho;
	pi = M_PARA/rho;
	x = -pi + (theta/2.0) + 0.5667;		/* vertical displacement, degrees */
	psi = sin(DEGTORAD(latitude))/cos(DEGTORAD(decs));
	psi = acos(psi);
	y = RADTODEG(asin(sin(DEGTORAD(x))/sin(psi)));	/* hor. displ, degrees */
	dt = (240.0*y)/cos(DEGTORAD(decs));		/* rise/set correction, sec */
	gsts += dt/3600.0;						/* since GST in hours */

	gst2ut(jd_date, gstr, &utr);
	gst2ut(jd_date, gsts, &uts);

	/* these may not be the right rise/set - could be we need tomorrow's */
	*jdr = jd_date + utr/24.0 + KLUDGE/1440.0;
	*jds = jd_date + uts/24.0 + KLUDGE/1440.0;
}

	/* return the phase of the moon at the given JD in the 'phase'
	   parameter, in units where

	     0.0  means moon is not illuminated at all (new)
	     0.5                illuminated on half of its extent (quarter)
	     1.0                fully illuminated (full)

	*/
	
void 
moon_phase(jd, phase)
double jd, *phase;
{
	double ra, dec, d, f, sm, slam;

	/* calculate the solar and lunar longitudes */
	sun_mlam(jd, &sm, &slam);
	moon_radec(jd, &ra, &dec);
	d = llong - slam;
	f = 0.5*(1.0 - cos(DEGTORAD(d)));
	*phase = f;
}