Moon Cycles


Description:
These two functions, and their support functions, can be used to determine the current phase of the moon given a particular date.
 
Code:
'--------------------------------------------------------------------
' Begin: Module modMoonCycle
'
' The function MoonPhase returns one of the Moon_Types for the
'  date passed.  Useful for creation of calendars, or when
'  you need to know the cycle of the moon.
'
' MoonPhaseDetail returns a double from 0 to .99 representing
'  the phase of the moon for a given day, where 0->.99 is a the
'  new moon through to the last quater
'
' The remaining functions are support functions, and should
'  not be called directly
'
'--------------------------------------------------------------------
' Adopted from pcal, which bares this copyright notice:
'
'  Routines to accurately calculate the phase of the moon
'
'  Originally adapted from "moontool.c" by John Walker, Release 2.0.
'
'   This routine (calc_phase) and its support routines were adapted
'   from phase.c (v 1.2 88/08/26 22:29:42 jef) in the program
'   "xphoon" (v 1.9 88/08/26 22:29:47 jef) by Jef Poskanzer and
'   Craig Leres. The necessary notice follows...
'
'   Copyright (C) 1988 by Jef Poskanzer and Craig Leres.
'
'   Permission to use, copy, modify, and distribute this software
'   and its documentation for any purpose and without fee is hereby
'   granted, provided that the above copyright notice appear in all
'   copies and that both that copyright notice and this permission
'   notice appear in supporting documentation.  This software is
'   provided "as is" without express or implied warranty.
'
'   These were added to "pcal" by RLD on 19-MAR-1991
'--------------------------------------------------------------------

Option Explicit

Public Enum Moon_Types
    Moon_Other = -1
    Moon_NewMoon = 0
    Moon_1stQuarter = 1
    Moon_FullMoon = 2
    Moon_3rdQuarter = 3
End Enum

'Astronomical constants.
Private Const epoch = 2444238.5 '1980 January 0.0

'Constants defining the Sun's apparent orbit.
Private Const elonge = 278.83354 'ecliptic longitude of the Sun
                                 'at epoch 1980.0
Private Const elongp = 282.596403 'ecliptic longitude of the Sun
                                  'at perigee
Private Const eccent = 0.016718 'eccentricity of Earth's orbit

'Elements of the Moon's orbit, epoch 1980.0.
Private Const mmlong = 64.975464 'moon's mean lonigitude at the epoch
Private Const mmlongp = 349.383063 'mean longitude of the perigee
                                   'at the epoch
Private Const mlnode = 151.950429 'mean longitude of the node at
                                  'the epoch

Private Const PI = 3.14159265358979
Private Const EPSILON = 0.000001

'Returns the phase of the moon as a number from 0.00 to 0.99
Public Function MoonPhaseDetail(InDate As Date) As Double

    Dim dateTemp As Date
    dateTemp = Int(InDate)
    MoonPhaseDetail = calc_phase(DatePart("m", dateTemp), _
        DatePart("d", dateTemp), DatePart("yyyy", dateTemp))

End Function

'Returns the approximate phase of the moon, or Moon_Other
' if none apply for the selected day
Public Function MoonPhase(InDate As Date) As Moon_Types

    Dim nPrev As Double
    Dim nCurr As Double
    Dim nNext As Double
    
    Dim dateTemp As Date
    dateTemp = Int(InDate) - 1
    nPrev = calc_phase(DatePart("m", dateTemp), _
        DatePart("d", dateTemp), DatePart("yyyy", dateTemp))
    
    dateTemp = Int(InDate)
    nCurr = calc_phase(DatePart("m", dateTemp), _
        DatePart("d", dateTemp), DatePart("yyyy", dateTemp))
    
    dateTemp = Int(InDate) + 1
    nNext = calc_phase(DatePart("m", dateTemp), _
        DatePart("d", dateTemp), DatePart("yyyy", dateTemp))
    
    MoonPhase = is_quarter(nPrev, nCurr, nNext)

End Function

'Begin of support functions:
Private Function torad(d As Double) As Double 'deg->rad
    torad = ((d) * (PI / 180#))
End Function
Private Function FNITG(x As Double) As Double 'cos from deg
    FNITG = (Sgn(x) * Int(Abs(x)))
End Function
Private Function fixangle(a As Double) As Double 'fix angle
    fixangle = ((a) - 360# * (Int((a) / 360#)))
End Function
Private Function todeg(d As Double) As Double 'rad->deg
    todeg = ((d) * (180# / PI))
End Function

'   julday -- calculate the julian date from input month, day, year
'       N.B. - The Julian date is computed for noon UT.
'
'       Adopted from Peter Duffett-Smith's book `Astronomy With Your
'       Personal Computer' by Rick Dyson 18-MAR-1991
Private Function julday(month As Long, day As Long, _
    year As Long) As Double

    Dim mn1 As Long, yr1 As Long
    Dim a As Double, b As Double, c As Double, d As Double
    Dim djd As Double

    mn1 = month
    yr1 = year
    If yr1 < 0 Then
        yr1 = yr1 + 1
    End If
    If month < 3 Then
        mn1 = month + 12
        yr1 = yr1 - 1
    End If
    
    If (year < 1582) Or ((year = 1582) And (month < 10)) Or _
        ((year = 1582) And (month = 10) And (day < 15#)) Then
        b = 0
    Else
        a = Int(yr1 / 100#)
        b = 2 - a + Int(a / 4)
    End If
    
    If yr1 >= 0 Then
        c = Int(365.25 * yr1) - 694025
    Else
        c = FNITG((365.25 * yr1) - 0.75) - 694025
    End If
    
    d = Int(30.6001 * (mn1 + 1))
    djd = b + c + d + day + 2415020#
    
    julday = djd

End Function

'  kepler - solve the equation of Kepler
Private Function kepler(m As Double, ecc As Double) As Double

    Dim e As Double, delta As Double

    e = torad(m)
    m = e
    Do
        delta = e - ecc * Sin(e) - m
        e = e - (delta / (1 - ecc * Cos(e)))
    Loop While (Abs(delta) > EPSILON)
    
    kepler = e

End Function

'  calc_phase - calculate phase of moon as a fraction:
'
'  The argument is the time for which the phase is requested,
'  expressed as the month, day and year.  It returns the phase of
'  the moon (0.0 -> 0.99) with the ordering as New Moon,
'  First Quarter, Full Moon, and Last Quarter.
'
'  Converted from the subroutine phase.c used by "xphoon.c" (see
'      above disclaimer) into calc_phase() for use in "moonphas.c"
'  by Rick Dyson 18-MAR-1991
Private Function calc_phase(month As Long, inday As Long, _
    year As Long) As Double

    Dim day As Double, n As Double, m As Double, ec As Double
    Dim lambdasun As Double, ml As Double, mm As Double
    Dim mn As Double, ev As Double, ae As Double, a3 As Double
    Dim mmp As Double, mec As Double, a4 As Double, lp As Double
    Dim v As Double, lpp As Double, np As Double, y As Double
    Dim x As Double, moonage As Double, pdate As Double

    'need to convert month, day, year into a Julian pdate
    pdate = julday(month, inday, year)

    'Calculation of the Sun's position.
    day = pdate - epoch 'date within epoch
    n = fixangle((360 / 365.2422) * day) 'mean anomaly of the Sun
    m = fixangle(n + elonge - elongp) 'convert from perigee
                                      'co-ordinates to epoch 1980.0
    ec = kepler(m, eccent) 'solve equation of Kepler
    ec = Sqr((1 + eccent) / (1 - eccent)) * Tan(ec / 2)
    ec = 2 * todeg(Atn(ec)) 'true anomaly
    lambdasun = fixangle(ec + elongp) 'Sun's geocentric ecliptic
                                      'longitude
    
    'Calculation of the Moon's position.
    ml = fixangle(13.1763966 * day + mmlong) 'Moon's mean longitude.
    mm = fixangle(ml - 0.1114041 * day - mmlongp) 'Moon's mean
                                                  'anomaly.
    mn = fixangle(mlnode - 0.0529539 * day) 'Moon's ascending node
                                            'mean longitude.
    ev = 1.2739 * Sin(torad(2 * (ml - lambdasun) - mm)) 'Evection.
    ae = 0.1858 * Sin(torad(m)) 'Annual equation.
    a3 = 0.37 * Sin(torad(m)) 'Correction term.
    mmp = mm + ev - ae - a3 'Corrected anomaly.
    mec = 6.2886 * Sin(torad(mmp)) 'Correction for the equation
                                   'of the centre.
    a4 = 0.214 * Sin(torad(2 * mmp)) 'Another correction term.
    lp = ml + ev + mec - ae + a4 'Corrected longitude.
    v = 0.6583 * Sin(torad(2 * (lp - lambdasun))) 'Variation.
    lpp = lp + v 'True longitude.
    
    'Calculation of the phase of the Moon.
    moonage = lpp - lambdasun 'Age of the Moon in degrees.
    
    calc_phase = (fixangle(moonage) / 360#)
End Function

'is_quarter - if current phase of moon coincides with quarter
' moon, return non negative integer, otherwise, return -1
Private Function is_quarter(prev As Double, curr As Double, _
    mnext As Double) As Long

    Dim quarter As Long
    Dim phase As Double, diff As Double
    
    'adjust phases for 1 -> 0 wraparound
    If curr < prev Then
        curr = curr + 1
    End If
    If mnext < prev Then
        mnext = mnext + 1
    End If

    is_quarter = -1

    ' if a quarter moon has occurred between "prev" and "next",
    ' return TRUE if "curr" is closer to it than "prev" or
    ' "next" is
    For quarter = 1 To 4
        phase = quarter / 4#
        If prev < phase Then
            If mnext > phase Then
                diff = Abs(curr - phase)
                If diff < phase - prev Then
                    If diff < mnext - phase Then
                        is_quarter = quarter Mod 4
                        Exit Function
                    End If
                End If
            End If
        End If
    Next

End Function

' End: Module modMoonCycle
'--------------------------------------------------------------------
 
Sample Usage:
 
    Dim nDate As Date
    nDate = #1/1/2001#
    Do Until year(nDate) <> 2001
        Select Case MoonPhase(nDate)
            Case Moon_NewMoon
                Debug.Print Format(nDate, "mm/dd/yyyy");
                Debug.Print " [  ] New moon"
            Case Moon_1stQuarter
                Debug.Print Format(nDate, "mm/dd/yyyy");
                Debug.Print " [ *] First quarter"
            Case Moon_FullMoon
                Debug.Print Format(nDate, "mm/dd/yyyy");
                Debug.Print " [**] Full moon"
            Case Moon_3rdQuarter
                Debug.Print Format(nDate, "mm/dd/yyyy");
                Debug.Print " [* ] Third quarter"
        End Select
        nDate = DateAdd("d", 1, nDate)
    Loop