/* standalone moon phase calculation */
/* clang -lm moon.c -o moon && ./moon */

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

#define PI 3.1415926535897932384626433832795
#define RAD (PI / 180.0)
#define SMALL_FLOAT (1e-12)

typedef struct {
	int year, month, day;
	double hour;
} TimePlace;

void
JulianToDate(TimePlace *now, double jd)
{
	long jdi, b, c, d, e, g, g1;
	jd += 0.5;
	jdi = jd;
	if(jdi > 2299160) {
		long a = (jdi - 1867216.25) / 36524.25;
		b = jdi + 1 + a - a / 4;
	} else
		b = jdi;
	c = b + 1524;
	d = (c - 122.1) / 365.25;
	e = 365.25 * d;
	g = (c - e) / 30.6001;
	g1 = 30.6001 * g;
	now->day = c - e - g1;
	now->hour = (jd - jdi) * 24.0;
	if(g <= 13)
		now->month = g - 1;
	else
		now->month = g - 13;
	if(now->month > 2)
		now->year = d - 4716;
	else
		now->year = d - 4715;
}

double
Julian(int year, int month, double day)
{
	/* Returns the number of julian days for the specified day. */
	int a, b = 0, c, e;
	if(month < 3) {
		year--;
		month += 12;
	}
	if(year > 1582 || (year == 1582 && month > 10) ||
		(year == 1582 && month == 10 && day > 15)) {
		a = year / 100;
		b = 2 - a + a / 4;
	}
	c = 365.25 * year;
	e = 30.6001 * (month + 1);
	return b + c + e + day + 1720994.5;
}

double
sun_position(double j)
{
	double n, x, e, l, dl, v, m2;
	int i;
	n = 360 / 365.2422 * j;
	i = n / 360;
	n = n - i * 360.0;
	x = n - 3.762863;
	if(x < 0) x += 360;
	x *= RAD;
	e = x;
	do {
		dl = e - .016718 * sin(e) - x;
		e = e - dl / (1 - .016718 * cos(e));
	} while(fabs(dl) >= SMALL_FLOAT);
	v = 360 / PI * atan(1.01686011182 * tan(e / 2));
	l = v + 282.596403;
	i = l / 360;
	l = l - i * 360.0;
	return l;
}

double
moon_position(double j, double ls)
{
	double ms, l, mm, n, ev, sms, z, x, lm, bm, ae, ec, d, ds, as, dm;
	int i;
	/* ls = sun_position(j) */
	ms = 0.985647332099 * j - 3.762863;
	if(ms < 0) ms += 360.0;
	l = 13.176396 * j + 64.975464;
	i = l / 360;
	l = l - i * 360.0;
	if(l < 0) l += 360.0;
	mm = l - 0.1114041 * j - 349.383063;
	i = mm / 360;
	mm -= i * 360.0;
	n = 151.950429 - 0.0529539 * j;
	i = n / 360;
	n -= i * 360.0;
	ev = 1.2739 * sin((2 * (l - ls) - mm) * RAD);
	sms = sin(ms * RAD);
	ae = 0.1858 * sms;
	mm += ev - ae - 0.37 * sms;
	ec = 6.2886 * sin(mm * RAD);
	l += ev + ec - ae + 0.214 * sin(2 * mm * RAD);
	l = 0.6583 * sin(2 * (l - ls) * RAD) + l;
	return l;
}

double
moon_phase(int year, int month, int day, double hour, int *ip)
{
	double j = Julian(year, month, (double)day + hour / 24.0) - 2444238.5;
	double ls = sun_position(j);
	double lm = moon_position(j, ls);
	double t = lm - ls;
	if(t < 0) t += 360;
	*ip = (int)((t + 22.5) / 45) & 0x7;
	return (1.0 - cos((lm - ls) * RAD)) / 2;
}

static void
nextDay(int *y, int *m, int *d, double dd)
{
	TimePlace tp;
	double jd = Julian(*y, *m, (double)*d);
	jd += dd;
	JulianToDate(&tp, jd);
	*y = tp.year;
	*m = tp.month;
	*d = tp.day;
}

int
main(int argc, char **argv)
{
	int y, m, d, m0, h, i;
	int ymax, mmax, dmax, hmax;
	int ymin, mmin, dmin, hmin;
	int begun = 0;
	double step = 1;
	double pmax = 0;
	double pmin = 1;
	double plast = 0;
	/* setup time */
	time_t seconds = time(NULL);
	struct tm zt = {0};
	struct tm *t = localtime(&seconds);
	if(t == NULL)
		t = &zt;
	/* set defaults */
	y = t->tm_year + 1900;
	m = t->tm_mon + 1;
	/* user input */
	printf("year[%d]: ", y);
	fflush(stdout);
	scanf("%d", &y);
	printf("month[%d]: ", m);
	fflush(stdout);
	scanf("%d", &m);
	/* ready */
	d = 1;
	m0 = m;
	printf("\nDate       Time   Phase Segment\n");
	for(;;) {
		double p;
		int ip;
		for(h = 0; h < 24; h += step) {
			p = moon_phase(y, m, d, h, &ip);
			if(begun) {
				if(p > plast && p > pmax) {
					pmax = p;
					ymax = y;
					mmax = m;
					dmax = d;
					hmax = h;
				} else if(pmax) {
					printf("%04d/%02d/%02d %02d:00          (fullest)\n",
						ymax,
						mmax,
						dmax,
						hmax);
					pmax = 0;
				}
				if(p < plast && p < pmin) {
					pmin = p;
					ymin = y;
					mmin = m;
					dmin = d;
					hmin = h;
				} else if(pmin < 1) {
					printf("%04d/%02d/%02d %02d:00          (newest)\n",
						ymin,
						mmin,
						dmin,
						hmin);
					pmin = 1.0;
				}
				if(h == 16) {
					printf("%04d/%02d/%02d %02d:00 %5.1f%%   (%d)\n",
						y,
						m,
						d,
						h,
						floor(p * 1000 + 0.5) / 10,
						ip);
				}
			} else
				begun = 1;
			plast = p;
		}
		nextDay(&y, &m, &d, 1.0);
		if(m != m0) break;
	}
	return 0;
}