/* Program Name: Moon Phases

Author: Ben Jacobs 

Compilation: gcc -o moon moon.c -lm (NOTE: must be linked to the math 

library, I don't know how you manage this in a non-UNIX environment)

Description: Gets user input for date & time and gives the current phase 

of the moon. The "fullness" of the moon is stored in the double precision

float "F". 

Comments: I don't take an credit for the moon phase math. I found it in a 

program written in QBASIC by J. Stein Carter who found the math in 

"Practical Astronomy with Your Calculator", 3rd Ed by Duffett-Smith, Peter. 

1988 Cambridge Univ. Press. All I did was translate it into C.

I realize it uses a pretty excessive amount of double precision floats but

that's how I found the calculations and I didn't want to mess with them 

that much.

Bug report: dooberwah@crosswinds.net

*/ 

#include 

#include 

#define PI 3.141592654

int day_of_year(int year, int month, int day);

int fullness(int d, int *numdays, double *F, double *x);

int main(void)

    {

     int month, day, year, leap;

     int hour, minute, second;

     int datenum, numdays;

     double F, x;

     /* Gets date from user. Returns -1 if user enters bad date */

     printf("Input date mm/dd/yyyy > ");

     scanf("%d/%d/%d", &month, &day, &year);

         if ((month < 1 || month > 12) || (day < 1 || day > 365) || year <= 0) {

         printf("Invalid date.\n");

         return -1;

         }

         /* Checks to see if year is a leap year */

         leap = year%4 == 0 && year%100 != 0 || year%400 == 0;

         printf("Input time in 24 hour formatt hh:mm:ss > ");

         scanf("%d:%d:%d", &hour, &minute, &second);

             if ((hour < 1 || hour > 24) || (minute < 1 || minute > 60) || (second < 1 || second > 60)) {

             printf("Invalid time.\n");

             return -1;

             }

             datenum = day_of_year(year, month, day);

             fullness(datenum, &numdays, &F, &x);

             if (F < 0.003) 

             printf("The moon is currently a new moon.\n");

             else if (F < 0.47)

             printf("The moon is in a first quarter moon.\n");

             else if (F > 0.47 && F < 0.53)

             printf("The moon is currently a half moon.\n");

             else if (F > 0.997)

             printf("The moon is currently a full moon.\n");

             else if (F > 0.53)

             printf("The moon is currently in the third quarter.\n");

             return 0;

        }

        /* Here's a structure with the number of days in each month. for leap and 

        non-leap years. It's like the monthname array in the way that it uses

        a 0 in the first place (or "Illegal month" in monthname) so that January

        is the first month instead of the 0th. */

            static char daytab[2][13] = {

             {0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}, 

             {0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}

        };

        /* this function figures the day number out of 365 for any given day in any 

        month. */

        int day_of_year(int year, int month, int day) 

            {

             int i, leap;

             leap = year%4 == 0 && year%100 != 0 || year%400 == 0;

             for (i = 1; i < month; i++)

             day += daytab[leap][i];

             return day;

        }

        int fullness(int d, int *numdays, double *F, double *x)

            {

             double nCounter, N, Mo, Ec, lamda, l, Mm, N65, Ev, Ae, A3, Mmp, A4, lp, V, lpp, D67;

             *numdays = d + nCounter / 2;

             N = *numdays * (360 / 365.242191);

             N = abs(N) - abs(N / 360) * 360; 

             Mo = N + 279.403303 - 282.768422; 

             if (Mo < 0)

             Mo = Mo + 360;

             Ec = 360 * .016713 * sin(Mo * PI / 180) / PI;

             lamda = N + Ec + 279.403303;

             if (lamda > 360)

             lamda = lamda - 360;

             l = 13.1763966 * *numdays + 318.351648;

             l = abs(l) - abs(l / 360) * 360;

             Mm = l - .1114041 * d - 36.34041;

             Mm = abs(Mm) - abs(Mm / 360) * 360;

             N65 = 318.510107 - .0529539 * *numdays;

             N65 = abs(N65) - abs(N65 / 360) * 360;

             Ev = 1.2739 * sin((2 * (1 - lamda) - Mm) * PI / 180);

             Ae = .1858 * sin(Mo * PI / 180);

             A3 = .37 * sin(Mo * PI / 180);

             Mmp = Mm + Ev - Ae - A3;

             Ec = 6.2886 * sin(Mmp * PI / 180);

             A4 = .214 * sin((2 * Mmp) * PI / 180);

             lp = l + Ev + Ec - Ae + A4;

             V = .6583 * sin((2 * (lp - lamda)) * PI / 180);

             lpp = lp + V;

             D67 = lpp - lamda;

             *F = .5 * (1 - cos(D67 * PI / 180));

             *x = (sin(D67 * PI / 180));

             return 0;

        }