/*******************************************************************************
** moonphase.c                                                                **
**                                                                            **
**     Phase of the Moon. Calculates the current phase of the moon.           **
**     Based on routines from `Practical Astronomy with Your Calculator',     **
**        by Duffett-Smith.                                                   **
**     Comments give the section from the book that particular piece          **
**        of code was adapted from.                                           **
**                                                                            **
**     -- Keith E. Brandt VIII 1984 (Obviously for Aztec C, as)               **
**                                                                            **
**     -- Conversion to SAS/C 5.1b, Clean Up and German Version               **
**        Correct usage of ANSI-C time routines, especially `gmtime'.         **
**        Andreas Scherer VII 1992                                            **
**                                                                            **
*******************************************************************************/

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

#define EPOCH    1983
#define EPSILONg  279.103035   /* solar ecliptic long at EPOCH            */
#define RHOg      282.648015   /* solar ecliptic long of perigee at EPOCH */
#define e           0.01671626 /* solar orbit eccentricity                */
#define lzero     106.306091   /* lunar mean long at EPOCH                */
#define Pzero     111.481526   /* lunar mean long of perigee at EPOCH     */
#define Nzero      93.913033   /* lunar mean long of node at EPOCH        */

#define RAD 0.01745329251994329651 /* PI / 180 */

char *_TZ = "CET-01CDT";

void main()
   {
   double days;      /* days since EPOCH */
   double phase;     /* percent of lunar surface illuminated */
   double phase2;    /* percent of lunar surface illuminated one day later */
   long   lo;        /* used in time call */
   int    i = EPOCH;
   int    i_phase;
   struct tm *pt;    /* ptr to time structure */
   double potm();
   int    ly();

#ifdef AMIGA
   tzset();
#endif /* AMIGA */

   time(&lo);        /* get system time */
   pt = gmtime(&lo); /* get ptr to gmt time struct */

/*******************************************************************************
** Calculate days since EPOCH.                                                **
*******************************************************************************/
   days = (pt->tm_yday + 1.) +
      (((double)pt->tm_hour + (pt->tm_min / 60.) +
       (pt->tm_sec / 3600.)) / 24.);
   while(i < pt->tm_year + 1900)
      days += 365. + (double)ly(i++);

   phase = potm(days);

#ifdef EUTSCH
   printf("\nDer Mond ist ");
#else
   printf("\nThe Moon is ");
#endif

   i_phase = (int)(phase + 0.5);

   if(i_phase == 100)

#ifdef EUTSCH
      printf("voll\n\n");
#else
      printf("Full\n\n");
#endif

   else {
      if(i_phase == 0)

#ifdef EUTSCH
         printf("neu\n\n");
#else
         printf("New\n\n");
#endif

      else {
         if(i_phase == 50)  {
            phase2 = potm(++days);
            if(phase2 > phase)

#ifdef EUTSCH
               printf("im ersten Viertel\n\n");
#else
               printf("at the First Quarter\n\n");
#endif

            else 

#ifdef EUTSCH
               printf("im letzten Viertel\n\n");
#else
               printf("at the Last Quarter\n\n");
#endif

            }
         else {
            if(i_phase > 50) {
               phase2 = potm(++days);
               if(phase2 > phase)

#ifdef EUTSCH
                 printf("zunehmend");
#else
                 printf("Waxing");
#endif

               else 

#ifdef EUTSCH
                  printf("abnehmend");
#else
                  printf("Waning");
#endif

#ifdef EUTSCH
               printf(", Halbmond (%1.0f%% voll)\n\n", phase);
#else
               printf(", Gibbous (%1.0f%% of Full)\n\n", phase);
#endif

               }
            else {
               if(i_phase < 50) {
                  phase2 = potm(++days);
                  if(phase2 > phase)

#ifdef EUTSCH
                     printf("zunehmend");
#else
                     printf("Waxing");
#endif

                  else

#ifdef EUTSCH
                     printf("abnehmend");
#else
                     printf("Waning");
#endif

#ifdef EUTSCH
                  printf(", Halbmond (%1.0f%% voll)\n\n", phase);
#else
                  printf(", Crescent (%1.0f%% of Full)\n\n", phase);
#endif

                  }
               }
            }
         }
      }
   }

/*******************************************************************************
** Calculate the Position Of The Moon to a given day fraction.                **
*******************************************************************************/
double potm(days)
double days;
   {
   double N,Msol,LambdaSol,Ev,Mmprime,lprime;
   double dtor(),range();

   N = range(360. * days / 365.2422) + EPSILONg;                /* sec 42 #03 */

   Msol = sin(dtor(range(N - RHOg)));                           /* sec 42 #04 */

   LambdaSol = range(N + 360./PI * e * Msol);                   /* sec 42 #05 */
                                                                /* sec 42 #06 */

   lprime = range(13.1763966 * days + lzero);                   /* sec 61 #04 */

   Mmprime = range(lprime - (0.1114041 * days) - Pzero);        /* sec 61 #05 */

   Ev = 1.2739 * sin(dtor(2. * (lprime - LambdaSol) - Mmprime)) /* sec 61 #07 */
      - 0.1858 * Msol;                                          /* sec 61 #08 */

   Mmprime += Ev - 0.3700 * Msol;                               /* sec 61 #09 */

   lprime += Ev                                                 /* sec 61 #10 */
      + 6.2886 * sin(dtor(Mmprime))                             /* sec 61 #11 */
      + 0.2140 * sin(dtor(2. * Mmprime));                       /* sec 61 #12 */

   return(50. * (1. - cos(dtor(lprime                           /* sec 61 #13 */
      + 0.6583 * sin(dtor(2. * (lprime - LambdaSol)))           /* sec 61 #14 */
      - LambdaSol))));                                          /* sec 63 #02 */
                                                                /* sec 63 #03 */
   }

/*******************************************************************************
** Returns 1 if LeapYear, 0 otherwise, according to Gregorian calendar.       **
*******************************************************************************/
int ly(yr)
int yr;
   {
   return(yr % 4 == 0 && yr % 100 != 0 || yr % 400 == 0);
   }

/*******************************************************************************
** Convert Degrees TO Radians                                                 **
*******************************************************************************/
double dtor(deg)
double deg;
   {
   return(deg * RAD);
   }