From:	MERC::"uunet!CRVAX.SRI.COM!RELAY-INFO-VAX" 23-JUN-1992 04:43:34.36
To:	info-vax@sri.com
CC:	
Subj:	Phases of the moon...Got it I think!

Info-Peoples,

     I think i got it!!!!  Thanx to Peter Quodling's sample code.  I converted
his DCL example to pascal.  I thought I would enter this as a summary.  I have
included the DCL and PASCAL versions.  I'm not sure if my program is right.
The only way i know to check it is to look out my window and check and then run
my program and compare.  Could some of you please extract this and check the
accuracy.  I would hate to put this in my program and it be wrong.  That
wouldn't be good.  Any Reply's would be greatly appreciated...  Here are
the programs...

                                      David Barbee
                                      STUBRBEE@EKU


$!  Updated for VMS V4.6, 16-Aug-1987
$! from Martin Minow's TECO program, which was from Knuth's "The
$! Art of Computer Programming", vol. 1.  Also, CACM, Apr 63.
$!
$  now      =   f$time()
$  month    =   f$cvtime(now,,"month")
$  day      =   f$cvtime(now,,"day")
$  daynumber    =   30*month-31 + day    ! Always <= real day-number (1-Jan =1)
$  date     =   f$cvtime(now,,"date")
$  basedate =   f$cvtime("1-JAN--1:0","absolute","date")
$find_day:
$  if f$cvtime(basedate+"+''DAYNUMBER'-",,"date") .eqs. date then goto got_date
$  daynumber    =   daynumber + 1
$  goto find_day
$got_date:
$  year     =   f$cvtime(now,,"year")
$  century  =   year/100 + 1
$  golden   =   year - ((year/19)*19) + 1
$  epact    =   golden*11+20 + ((8*century+5)/25-5) - (3*century/4-12)
$  epact    =   epact - ((epact/30)*30)
$  if epact .eq. 25 .and. golden .gt. 11 then epact = 26
$  if epact .eq. 24 then epact = 25
$! Calculate the phase of the moon.
$!  The moon's period is 29.5 days; we need to /8, and then /2,
$!  and the period already has a divisor of 2, so we *64, giving
$!  the moon's period expressed in sixty-fourths of a day.
$  multiplier   =   64
$  period   =   29*multiplier+multiplier/2      ! =29.5 * 64
$  eighth   =   period/8       ! =29.5 * 8      ! duration of one phase
$  sixteenth    =   eighth/2    ! =29.5 * 4      ! used for rounding
$  phase    =   (daynumber+epact-1)*multiplier+sixteenth
$  phase    =   ((phase - ((phase/period)*period))/eighth) .and. %O7
$  phase_names  =   "new,waxing crescent,in the first quarter,"+ -
            "waxing gibbous,full,waning gibbous,"+ -
                    "in the last quarter,waning crescent"
$ write sys$output "The moon is ",f$element(phase,",",phase_names),"."


Program Moon_Phase(Input, Output);

Const
  Multiplier = 64;

Var
  Date        : Timestamp;
  Day_Of_Year,
  Century,
  Golden,
  Epact,
  Period,
  Eighth,
  Sixteenth,
  Phase       : Integer;
  Phase_Names : Array [0..7] Of Varying [20] of Char;

Value
  Phase_Names := ('New', 'Waxing Crescent', 'First Quarter', 'Waxing Gibbous',
                  'Full','Waning Gibbous', 'Last Quarter', 'Waning Crescent');

[Asynchronous]
Function Syb$Day_Of_Week(Day,
                         Month,
                         Year  : Integer := 0) : Integer;

Var
  m       : integer;                                        (* Adjusted month, l
   et year start at march 1, so leap day only at end *)
  y       : integer;                                        (* Year in current c
   entury, adjusted -1 if month was jan or feb       *)
  c       : integer;                                        (* Century number
                                                     *)
  Default : TimeStamp;

Begin
  If (((Year <= 1752) and (Month <= 2) and (Day < 2)) And (Year <> 0))
    Then
      Syb$Day_Of_Week := 0
    Else
      Begin
        If Day = 0
          Then
            Begin
              GetTimeStamp(Default);
              Day   := Default.Day;
              Month := Default.Month;
              Year  := Default.Year;
            End;
        y := year;
        m := month - 2;
        if m < 1 then
          Begin
            m := m + 12;
            y := pred(y);
          end;
        c := y div 100;
        y := y mod 100;
        Syb$Day_Of_Week := succ(((260 * m - 15) div 100 + day + y + y div 4 + c
   div 4 - 2 * c) mod 7);
      end;
End;

[Asynchronous]
Function Syb$Day_Of_Year(Day,
                         Month,
                         Year  : Integer := 0) : Integer;

Var
  Days_In_Month    : [Static] Array [1..12] of Integer;
  Default          : TimeStamp;
  Loop             : Integer;
  Temp             : Integer := 0;

Value
  Days_In_Month := (31,28,31,30,31,30,31,31,30,31,30,31);

Begin
  If Day = 0
    Then
      Begin
        GetTimeStamp(Default);
        Day   := Default.Day;
        Month := Default.Month;
        Year  := Default.Year;
      End;
  If (Year Mod 4) = 0 then Days_In_Month[2] := 29;
  For Loop := 1 to Month - 1 Do
    Temp := Temp + Days_In_Month[Loop];
  Temp := Temp + Day;
  Syb$Day_Of_Year := Temp;
End;

Begin
  GetTimeStamp(Date);
  Day_Of_Year := Syb$Day_Of_Year(Date.Day,Date.Month,Date.Year);
  Century     := Date.Year Div 100 + 1;
  Golden      := Date.Year - ((Date.Year Div 19) * 19) + 1;
  Epact       := Golden * 11 + 20 + ((8 * Century + 5) Div 25 - 5) - (3 * Centur
   y Div 4 - 12);
  Epact       := epact - ((epact Div 30) * 30);
  If ((Epact = 25) And (Golden > 11)) Then Epact := 26;
  If (Epact = 24) Then epact := 25;
  Period      := 29 * Multiplier + Multiplier Div 2;
  Eighth      := Period Div 8;
  Sixteenth   := eighth Div 2;
  Phase       := (Day_Of_Year + Epact - 1) * Multiplier + Sixteenth;
  Phase       := (( Phase - ((Phase Div Period) * Period)) Div Eighth);
  Repeat
    If Phase > 7 then Phase := Phase - 8;
  Until Phase <= 7;
  Writeln(Phase_Names[Phase]);
End.

