[comp.lang.c] Calculating Phases of the Moon

hicks@walker-emh.arpa (Gregory Hicks COMFLEACTS) (02/14/88)

While I was an under-grad at the University of Utah many moons ago, I
used a program called 'calendar' that allowed you to make a calendar that
started on whatever day of the week you wanted your week to start on (Mon,
Tue, etc), allowed you to list important days by specifying the date, or day
of the week and frequency (ie: every other Monday or Every Monday).

This calendar program figured out what day of the week the year started on
by calculating what date in March/April Easter fell on (this is based on the
phase of the moon...) and then counting backwards an appropriate number of
days (calculating back really...)

This calculation was based on something called 'The Golden Number' and was
part of a programming problem given to beginning programming students...

The algorithm was correct after the calendar changed from the Julian
Calendar to the Gregorian Calendar... (Before that, you had to take into
account the number of days that were 'lost' in the change...

I would like very much to get my hands on this program because...  First I
have a requirement to calculate phases of the moon, second, I
want to port it to another system, and LASTLY, I need the verbage that
describes the 'History' of the Gregorian Calendar...

Anybody have any idea where it might be found?  The guy that wrote it
(I THINK) was Carl M. Ellison.  Unfortunately, I don't know where to fing
him anymore because it HAS BEEN over 17 years since I left there...

Regards,
Gregory Hicks
<hicks@walker-emh.arpa>

kenny@uiucdcsb.cs.uiuc.edu (02/16/88)

For the Lilius-Clavius algorithm for the date of Easter in the
Gregorian calendar, see Knuth, Donald E. _The Art of Computer
Programming:_ Volume 1: Fundamental Algorithms.  2d ed., Reading,
Massachusetts, Addison-Wesley, 1973, pp. 155-156. DON'T FORGET
EXERCISE 15!

Other references:

Communications of the ACM 5 (1962) 209-210, gives an algorithm which
includes the Julian calendar. and discusses the algorithm given by
Victor of Aquitaine.   Also look at O'Beirne, _Puzzles and Paradoxes._
London: Oxford University Press, 1965, chapter 10.

For the phase of the astronomical moon, rather than the full moon, try
the following program:
------------------------------------------------------------------------
/****************************************************************************
 pom.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

! modified by awpaeth@watcgl, December 1987 for sysV compatability
****************************************************************************/

#include <stdio.h>
#ifndef SYSV
#include <sys/time.h>
#else
#include <time.h>
#endif
#include <math.h>
#define PI         3.141592654
#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 */

main()  {

double dtor();
double adj360();
double potm();

long *lo = (long *) calloc (1, sizeof(long)); /* used by time calls */
struct tm *pt; /* ptr to time structure */

double days;   /* days since EPOCH */
double phase;  /* percent of lunar surface illuminated */
double phase2; /* percent of lunar surface illuminated one day later */
int i = EPOCH;

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

/* calculate days since EPOCH */
days = (pt->tm_yday +1) + ((pt->tm_hour + (pt->tm_min / 60.0)
       + (pt->tm_sec / 3600.0)) / 24.0);
while (i < pt->tm_year + 1900)
   days = days + 365 + ly(i++);

phase = potm(days);
printf("The Moon is ");
if ((int)(phase + .5) == 100) {
   printf("Full\n");
   }
else if ((int)(phase + 0.5) == 0) 
   printf("New\n");
else if ((int)(phase + 0.5) == 50)  {
   phase2 = potm(++days);
   if (phase2 > phase)
      printf("at the First Quarter\n");
   else 
      printf("at the Last Quarter\n");
   }
else if ((int)(phase + 0.5) > 50) {
   phase2 = potm(++days);
   if (phase2 > phase)
      printf("Waxing ");
   else 
      printf("Waning ");
   printf("Gibbous (%1.0f%% of Full)\n", phase);
   }
else if ((int)(phase + 0.5) < 50) {
   phase2 = potm(++days);
   if (phase2 > phase)
      printf("Waxing ");
   else
      printf("Waning ");
   printf("Crescent (%1.0f%% of Full)\n", phase);
   }
}

double potm(days)
double days;
{
double N;
double Msol;
double Ec;
double LambdaSol;
double l;
double Mm;
double Ev;
double Ac;
double A3;
double Mmprime;
double A4;
double lprime;
double V;
double ldprime;
double D;
double Nm;

N = 360 * days / 365.2422;  /* sec 42 #3 */
adj360(&N);

Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
adj360(&Msol);

Ec = 360 / PI * e * sin(dtor(Msol)); /* sec 42 #5 */

LambdaSol = N + Ec + EPSILONg;       /* sec 42 #6 */
adj360(&LambdaSol);

l = 13.1763966 * days + lzero;       /* sec 61 #4 */
adj360(&l);

Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */
adj360(&Mm);

Nm = Nzero - (0.0529539 * days);     /* sec 61 #6 */
adj360(&Nm);

Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */

Ac = 0.1858 * sin(dtor(Msol));       /* sec 61 #8 */
A3 = 0.37 * sin(dtor(Msol));

Mmprime = Mm + Ev - Ac - A3;         /* sec 61 #9 */

Ec = 6.2886 * sin(dtor(Mmprime));    /* sec 61 #10 */

A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */

lprime = l + Ev + Ec - Ac + A4;      /* sec 61 #12 */

V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */

ldprime = lprime + V;                /* sec 61 #14 */

D = ldprime - LambdaSol;             /* sec 63 #2 */

return (50 * (1 - cos(dtor(D))));    /* sec 63 #3 */
}

ly(yr)
int yr;
{
/* returns 1 if leapyear, 0 otherwise */
return (yr % 4 == 0 && yr % 100 != 0 || yr % 400 == 0);
}

double dtor(deg)
double deg;
{
/* convert degrees to radians */
return (deg * PI / 180);
}

double adj360(deg)
double *deg;
{
/* adjust value so 0 <= deg <= 360 */
do if (*deg < 0)
   *deg += 360;
else if (*deg > 360)
   *deg -= 360;
while (*deg < 0 || *deg > 360);
}
------------------------------------------------------------------------
Kevin Kenny			UUCP: {ihnp4,pur-ee,convex}!uiucdcs!kenny
Department of Computer Science	ARPA: kenny@B.CS.UIUC.EDU (kenny@UIUC.ARPA)
University of Illinois		CSNET: kenny@UIUC.CSNET
1304 W. Springfield Ave.
Urbana, Illinois, 61801		Voice: (217) 333-8740