[net.astro] Phase of Moon from Unix time

doug@escher.UUCP (Douglas J Freyburger) (02/24/86)

I am looking for a routine that will return the phase of
the moon (in floating, a 0-28 integer, or  whatever) from a
Unix binary time gotten from 'long time()' or any of the
other time routines.  I am only interested in the level of
accuracy given on calendars with the little pictures of the
new moon, first quarter, full moon, and third quarter on
them, so timezone and daylight savings corrections are
optional.  Alternately, anyone who has already built a
phase of the moon feature into your 'calendar' utility,
would you be willing to share the code?
	(-)NX in advance,
	Doug Freyburger
-- 

Doug Freyburger		DOUG@JPL-VLSI,
JPL Mail Stop 23	escher!doug, escher!teleop!doug
Pasadena, CA 91109	etc.

<Generic Disclaimer>

paul@zoom.UUCP (20 Proprietor) (02/26/86)

A program came across the net some time in the recent past.  It was written
by "ArchMach", converted from PL/I by Brian Hess and "extensively cleaned
up" by Rich $alz.  It's called `pom' and I'm sending it to you under separate
cover.  If there's enough interest, you may want to post it to the world.  It
is certainly accurate enough for most, and can be quickly modified to give
you a floating point number representing the phase (rather than a character
string).  I have a multitude of (sometimes ridiculously accurate) programs
for the determination of the ephemerides.  Such topics as:

	* Time conversion routines - UN*X time, UT, ET, UT adjusted for
		longitude, Sidereal Time (apparent, mean), etc., all
		with routines to facilitate cross-conversion.
	* Nutation calculations - based on 1984 IAU Theory of Nutation
		(yes, I DID enter all those coefficients!)
	* Sun calculations
	* Moon calculations
	* Solution of Kepler's equation
	* Comet ephemeris routines (including special programs for the
		oldest sensation...comet Halley)

By the way, they're all written in `C' and run on my 3b2/300 (without a
floating point coprocessor, they actually jog :').  Anyone interested can
send me mail and I'd be glad to further explain, etc.  Given the sheer volume
(about 4k lines), I'd rather not post to the public.  At any rate, feel free
to reply if you have any interest in astronomical calculations.

Who says JPL's cornered the market on space computations :'> ?

								Paul Ruel
							ihnp4!ariel!zoom!paul

graham@convex.UUCP (02/27/86)

> I am looking for a routine that will return the phase of
> the moon (in floating, a 0-28 integer, or  whatever) ...

I'd like to have this code too.

Marv Graham; Convex Computer Corp. {allegra,ihnp4,uiucdcs,ctvax}!convex!graham

dar@telesoft.UUCP (David Reisner @favorite) (03/06/86)

> > I am looking for a routine that will return the phase of
> > the moon (in floating, a 0-28 integer, or  whatever) ...

This was in net.sources sometime in the past.  It had some bugs in it,
which I have fixed (someone's C compiler does not perform the correct
implicit conversions in expression evaluation).  It does a little more
than just return the phase...

-David
sdcsvax!telesoft!dar


/*  The phase of the moon, for your safety and convenience.
**  (compile cc pom.c -lm	hjs)
**
**  Stolen from ArchMach & converted from PL/I by Brian Hess.
**  Extensively cleaned up by Rich $alz.
**
**  If you can figure out how this program works, then YOU deserve
**  to be working on it, sucker!  Here's a hint:  The epoch is 13:23,
**  10/1/78.
*/
/*  Change Log
    25jun85	dar	fixed day/hour comma printing; calculate Fraction w/dbls
*/


#include <math.h>
#include <sys/types.h>
#include <time.h>


/* Globals. */
long		 Day;
long		 Hour;
long		 Minute;
double		 Fraction;	/*-dar 0.0=full, 0.25=last qrtr, 0.5=new, */
				/*-dar 0.75=first qrtr, 1.0=full */
char		*Moon[] =
{
    "new",
    "first quarter of the",
    "full",
    "last quarter of the"
};


/* Linked in later. */
char		*rindex();
time_t		 time();
struct tm	*localtime();


#define LINES		24
#define WIDTH		80
#define CENTER		((WIDTH - 2 * LINES) / 2)
#define BRIGHT		'@'
#define LEDGE		'('
#define REDGE		')'
#define FULL		0.5
#define TwoPi		(2 * 3.14159)
#define ZERO		0.03

#define Plural(X)	if (X != 1) printf("s");
/*!*/
long
Calculate()
{
    register struct tm	*tm;
    register long	 Length;
    register long	 Phase;
    register long	 DeltaH;
    register long	 Delta;
    register long	 offset;
    time_t		 tick;
    long		 julian;
    long		 year;
    long		 hour;
    long		 minute;

    /*-dar
      Assume that the following code computes the current year, hour, and
      minute.
    */
    tick	= time((time_t *)0);
    tm		= localtime(&tick);
    julian	= tm->tm_yday + 1;
    year	= tm->tm_year - 78;
    hour	= tm->tm_hour;
    minute	= tm->tm_min;

    /*-dar
      Assume that
		Length is the length of a lunar cycle (units?),
		0 <= Phase <= Length is some portion of the cycle Length
			and is computed by position in epoch (Delta) mod
			cycle Length.
    */
    Length	= (double)2551 / 60 * 1000 + (double)443 / 60;
    offset	= ((year * 365L + julian) * 24L + hour) * 60L + minute;
    Delta	= offset - (273L * 24L + 13L) * 60L + 23L;
    Phase	= Delta - (Delta / Length) * Length;

    /*-dar
      Phase / Length should be between 0 and 1, and should be the desired
      Fraction of the moon to display.

      C defines <float> = <long> / <long> as an integer division, a
      conversion (of the result of the division) to float, and an
      assignment.  We need float division to get a non-zero result, as
      desired.
    */
    Fraction	= (double)Phase / (double)Length;

    Length	= Length / 4;
    Phase	= Phase / Length;
    Delta	= Delta - (Delta / Length) * Length;
    DeltaH	= Delta / 60;
    Minute	= Delta - DeltaH * 60;
    Day		= DeltaH / 24;
    Hour	= DeltaH - (Day * 24);
    /*printf("recomputed Fraction would be %lf\n", Phase / Length ); <<<< -dar*/
    return(Phase);
}
/*!*/
int
CharPos(x)
    double		x;
{
    register int	i;

    i = x * LINES + 0.5;
    if ((i += LINES + CENTER) < 1)
	i = 1;
    return(i);
}


void
Draw()
{
    register char	*p;
    register int	 i;
    register int	 end;
    register double	 y;
    register double	 cht;
    register double	 squisher;
    register double	 horizon;
    register double	 terminator;
    char		 Buffer[256];

    /* Clear screen? */

    if (Fraction < FULL)
	squisher = cos(TwoPi * Fraction);
    else
	squisher = cos(TwoPi * (Fraction - FULL));

    cht = (double)2.0 / (LINES - 6.0);
    for (y = 0.93; y > -1.0; y -= cht)
    {
	for (i = sizeof Buffer, p = Buffer; --i >= 0; )
	    *p++ = ' ';
	horizon = sqrt(1.0 - y * y);
	Buffer[    CharPos(-horizon)]	= LEDGE;
	Buffer[i = CharPos( horizon)]	= REDGE;
	Buffer[++i]			= '\0';
	terminator = horizon * squisher;
	if (Fraction > ZERO && Fraction < (1.0 - ZERO))
	{
	    if (Fraction < FULL)
	    {
		i   = CharPos( terminator);
		end = CharPos( horizon);
	    }
	    else
	    {
		i   = CharPos(-horizon);
		end = CharPos( terminator);
	    }
	    while (i <= end)
		Buffer[i++] = BRIGHT;
	}
	printf(" %s\n", Buffer);
    }
}
/*!*/
main(ac)
    int			 ac;
{
    register long	 Phase;

    Phase = Calculate();
    if (ac == 1)
	Draw();

    printf("\nIt is ");
    if (Day == 0 && Hour == 0 && Minute == 0)
	printf("exactly");
    else
    {
	if (Day)
	{
	    printf("%ld day", Day);
	    Plural(Day);
	    if (Hour | Minute)		/*-dar was Hour | Day */
		printf(", ");
	}
	if (Hour)
	{
	    printf("%ld hour", Hour);
	    Plural(Hour);
	    if (Minute)
		printf(", ");
	}
	if (Minute)
	{
	    printf("%ld minute", Minute);
	    Plural(Minute);
	}
	printf(" since");
    }
    printf(" the %s moon.\n", Moon[Phase]);

    exit(0);
}