VBScript to calculate moon age and phase

Windows Scripting Host

This subroutine calculates the age and phase of the moon for the supplied date string.

Sample output:

Moon's age (days): 2
Moon's phase: waxing crescent
Distance (Earth radii): 63.47
Ecliptic latitude (degrees): -4.31
Ecliptic longitude (degrees): 21.27

Notes: the description of the moon's phase is based on commonly
used terms and is not represented as accurate by any other
standard. It is a rough approximation.

The formulas used herein were found on the WWW. Unfortunately
I no longer have the original citation. My thanks to the person who
made them available on the www.

'Code begins here

' Moonphase.vbs
' 3/16/2002
' Copyright Paul R. Sadowski <aa089(at)bfn.org>
'===============================================
option explicit
MoonAge(Date)

Sub MoonAge(strDate)
dim tmp, arrtmp, strmonth, strday, stryear, strsep
dim yy, y, m, d, p2, mm, k1, k2, k3, l, j, v, ip, ag, dp, di, la
dim lo, po, rp, np, phase, thisphase

Tmp = InStr(1, strDate, "/")
if Tmp <> 0 then
strSep = "/"
else
Tmp = InStr(1, strDate, "-")
strSep = "-"
end if

arrTmp = Split(strDate, strSep,3)
strMonth = arrTmp(0)
strDay = arrTmp(1)
strYear = arrTmp(2)

Y = strYear
M = strMonth
D = strDay
P2=2*3.14159

YY=Y-int((12-M)/10)
MM=M+9
if MM>=12 then
MM=MM-12
end if
K1=int(365.25*(YY+4712))
K2=int(30.6*MM+.5)
K3=int(int((YY/100)+49)*.75)-38

' JD for dates in Julian calendar
J=K1+K2+D+59
if J>2299160 then
' For Gregorian calendar
J=J-K3
end if

' J is the Julian date at 12h UT on day in question

' Calculate illumination (synodic) phase
V=(J-2451550.1)/29.530588853
V=V-int(V)
if V<0 then
V=V+1
end if
IP=V

' Moon's age in days
AG=IP*29.53

' Convert phase to radians
IP=IP*P2

' Calculate distance from anomalistic phase
V=(J-2451562.2)/27.55454988
V=V-int(V)
if V<0 then
V=V+1
end if
DP=V
DP=DP*P2: ' Convert to radians
DI=60.4-3.3*COS(DP)-.6*COS(2*IP-DP)-.5*COS(2*IP)

' Calculate latitude from nodal (draconic) phase
V=(J-2451565.2)/27.212220817
V=V-int(V)
if V<0 then
V=V+1
end if
NP=V

' Convert to radians
NP=NP*P2
LA=5.1*SIN(NP)

' Calculate longitude from sidereal motion
V=(J-2451555.8)/27.321582241
' Normalize values to range 0 to 1
V=V-int(V)
if V<0 then
V=V+1
end if

RP=V
LO=360*RP+6.3*SIN(DP)+1.3*SIN(2*IP-DP)+.7*SIN(2*IP)

' phases from
'http://home.hiwaay.net/~krcool/Astro/moon/moonphase/
Phase = Array ( "new", "waxing crescent", "in its first quarter", _
"waxing gibbous", "full", "waning gibbous", _
"in its last quarter", "waning crescent" )
select case int(AG)
case 0, 29
ThisPhase = 0
case 1, 2, 3, 4, 5, 6
ThisPhase = 1
case 7
ThisPhase = 2
case 8, 9, 10, 11, 12, 13
ThisPhase = 3
case 14
ThisPhase = 4
case 15, 16, 17, 18, 19, 20, 21
ThisPhase = 5
case 22
ThisPhase = 6
case 23, 24, 25, 26, 27, 28
ThisPhase = 7
end select

wscript.echo "Moon's age (days): " & Int(AG)
wscript.echo "Moon's phase: " & Phase(ThisPhase)
wscript.echo "Distance (Earth radii): " & FormatNumber(DI,2)
wscript.echo "Ecliptic latitude (degrees): " & FormatNumber(LA,2)
wscript.echo "Ecliptic longitude (degrees): " & FormatNumber(LO,2)

End Sub