/* This file is moonphase.c.  It is one of the three .c files that must be
compiled to create the executable program moon. */

#if DEBUG
#include <stdio.h>
#endif
#include <math.h>
#include "moon.h"

/* Synodic month is number of days from one new moon to the next. */
/* Anomalistic month is number of days from moon's perigee to the next. */
/* Anomalistic year is number of days from earth's perihelion to the next. */
#define	MSYNODIC	29.530589
#define	MANOMALISTIC	27.55455
#define	YANOMALISTIC	365.2596413105

/* DSIN(theta) is the sine of the angle theta, which is in degrees.
Multiply theta by PI/180 to convert it to radians. */
#define DSIN(theta)	sin((PI * (theta)) / 180.0)


/* Return the current phase of the moon in radians.
	0	== new moon
	pi/2	== first quarter
	pi	== full moon
	3*pi/4	== last quarter
*/

double phase(int year, int month, int day, int hour, int minute, int second)
{
	int gregorian;	/* 1 for gregorian, 0 for pre-gregorian */
	double julianday;		/* holds a large number of days */
	double tillnew;	/* number of days until nect new moon */
	double angle;
	double synodics;	/* number of synodic months */

#ifdef DEBUG
	printf ("year, month, day, hour, min, sec: %d %d %d %d %d %d\n",
		year, month, day, hour, minute, second);
#endif
	/* Julian calendar before October 5, 1582; Gregorian calendar after. */
	if (year != 1582) {
		gregorian = year > 1582;
	} else if (month != 10) {
		gregorian = month > 10;
	} else {
		gregorian = day >= 5;
	}

	if (month <= 2) {
		--year;
		month += 12;
	}

	/* Julian days count from January 1, 4713 BC. */
	julianday = 365.0 * 4713.0 + 750.0 +
		(int)(365.25 * year) +
		(int)(30.6001 * (month + 1)) +
		day +
		hour / 24.0 +
		minute / 24.0 / 60.0 +
		second / 24.0 / 60.0 / 60.0;

	if (gregorian) {
		int centuries = year / 100;
		julianday -= centuries;
		julianday += centuries/4;
		julianday += 2;
	}
#ifdef DEBUG
	printf ("%.5f julianday\n", julianday);
#endif

	/* Number of synodic months elapsed from 4713 BC. */
	synodics = julianday / MSYNODIC + .67094;

	/* Number of synodic months until next new moon. */
	synodics -= (int)synodics;
	synodics = 1 - synodics;

	/* Number of days until next new moon. */
	tillnew = synodics * MSYNODIC;
#ifdef DEBUG
	printf ("tillnew == %.5f\n", tillnew);
#endif

	/* Julian day of next new moon. */
	julianday += tillnew;

	/* degrees in solar orbit at next new moon. */
	angle = (360.0 / YANOMALISTIC) * julianday;

	/* Number of days until next new moon, with fine correction. */
	tillnew += .1743 * DSIN(angle + 73.63);

	/* Convert Julian day of next new moon to degs in lunar orbit*/
	angle = (360.0 / MANOMALISTIC) * julianday + 271.5;

	/* Number of days until next new moon, with fine correction. */
	tillnew -= .4089 * DSIN(angle);

	/* Number of days until next new moon, with fine correction. */
	tillnew += .0161 * DSIN(2.0 * angle);
#ifdef DEBUG
	printf ("%.5f days till next new moon\n", tillnew);
	printf ("%.5f days since last new moon\n", MSYNODIC - tillnew);
#endif
	return 2 * PI * (MSYNODIC - tillnew) / MSYNODIC;
}