VBScript to calculate moon age and phase

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