/* standalone moon phase calculation */
#include <stdio.h>
#include <math.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;
long 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,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;
double 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;
double d;
double 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)
{
/*
Calculates more accurately than Moon_phase , the phase of the moon at
the given epoch.
returns the moon phase as a real number (0-1)
*/
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;
}
main()
{
int y, m, d;
int m0;
int h;
int i;
double step = 1;
int begun = 0;
double pmax = 0;
double pmin = 1;
int ymax, mmax, dmax, hmax;
int ymin, mmin, dmin, hmin;
double plast = 0;
printf("tabulation of the phase of the moon for one month\n\n");
printf("year: "); fflush(stdout);
scanf("%d", &y);
printf("month: "); 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;
}