*

1 files changed, 216 insertions, 0 deletions

diff --git a/c/moon.c b/c/moon.c
new file mode 100644
index 0000000..d564bd1
--- /dev/null
+++ b/c/moon.c
@@ -0,0 +1,216 @@
+/* standalone moon phase calculation */
+/* clang -lm moon.c -o moon && ./moon */
+/* via http://wiki.xxiivv.com/etc/moon.c.txt */
+
+#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;
+	time_t seconds = time(NULL);
+	struct tm zt = {0};
+	struct tm *t = localtime(&seconds);
+	if(t == NULL)
+		t = &zt;
+	printf("year[%d]: ", t->tm_year + 1900);
+	fflush(stdout);
+	scanf("%d", &y);
+	printf("month[%d]: ", t->tm_mon + 1);
+	fflush(stdout);
+	scanf("%d", &m);
+	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;
+}
\ No newline at end of file