/* standalone moon phase calculation */ /* clang -lm moon.c -o moon && ./moon */ #include #include #include #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; }