| '--------------------------------------------------------------------
' 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
'--------------------------------------------------------------------
|