VBScript to calculate sunrise and sunset

VBScript To Calculate Sunrise and Sunset

This VBscript is a port of the Perl Astro:Sunrise package by Ron Hill, rkhill(at)pacbell.net

The function SunTimes returns the sunrise and sunset times for any given date at a specific location.

Usage:
  SunTimes (strYear, strMonth, strDay, longitude, latitude, TZ, isdst)
Where:
  strYear is a string containing the year value of the date at question
  strMonth is a string containing the month of the date at question
  strDay is a string containing the day of the month of the date at question
  longitude is a numerical value representing the longitude of the local site (negative values for West)
  latitude is a numerical value representing the latitude of the local site
  TZ is a numerical offset from UTC at the local site
  isdst is a numerical value indicating whether Daylight Savings Time is in effect (1) or not (0).

SunTimes returns a comma seperated string indicating the sunrise and sunset times (24 hour clock), repectively. For example:
6:32,18:16

SunTimes is dependant on a number of support functions included below.

The following illustrates how to call the function to return the sunrise and sunset values for Buffalo, NY for the current date.

Const longitude = -78.79119894
Const latitude = 42.91730004
Const tzoff = -5

ThisTime = SunTimes(Year(Now), Month(Now), Day(Now), longitude, latitude, tzoff, 0)
Tmp = split(ThisTime, ",")
strOut = "Sunrise: " & Tmp(0) & vbCRLF & "Sunset: " & Tmp(1)
wscript.echo strOut

This code should adapt easily to either ASP or Visual Basic.

'Code begins here

'SunTimes.vbs
'=========================================================================
' Other Constants & Defs
Const pi = 3.1415926535897932384626433832795
RADEG=180/pi
DEGRAD=pi/180

'=========================================================================
' Main function
Function SunTimes (strYear, strMonth, strDay, longitude, latitude, TZ, isdst)
Dim d, n, i, w, m, l, e, e1, a, xv, yv, v, xs, ys, xe, ecl, lonsun, ye, ze, ra, dec, h
Dim GMST0, UT_Sun_in_south, LHA, hour_rise, hour_set, min_rise, min_set

'calculate days since 2000 jan 1
d = ( 367 * (strYear) - int((7 * ((strYear) + (((strMonth) + 9) / 12))) / 4) + int((275 * (strMonth)) / 9) + (strDay) - 730530)

' Orbital elements of the Sun:
N = 0.0
i = 0.0
w = 282.9404 + 4.70935E-5 * d

a = 1.000000
e = 0.016709 - 1.151E-9 * d

M = 356.0470 + 0.9856002585 * d
M = rev(M)

  ecl = 23.4393 - 3.563E-7 * d
  L = w + M
  if (L < 0 OR L > 360) then
    L = rev(L)
  end if

' position of the Sun
  E1 = M + e*(180/pi) * sind(M) * ( 1.0 + e * cosd(M) )
  xv = cosd(E1) - e
  yv = sqrt(1.0 - e * e) * sind(E1)

  v = atan2d(yv, xv)
  r = sqrt(xv * xv + yv * yv)
  lonsun = v + w
  if (lonsun < 0 OR lonsun > 360) then
    lonsun = rev(lonsun)
  end if
  xs = r * cosd(lonsun)
  ys = r * sind(lonsun)
  xe = xs
  ye = ys * cosd(ecl)
  ze = ys * sind(ecl)
  RA = atan2d(ye, xe)
  Dec = atan2d(ze, (sqrt((xe * xe)+(ye * ye))))
  h=-0.833

  GMST0 = L + 180
  if (GMST0 < 0 OR GMST0 > 360) then
    GMST0 = rev(GMST0)
  end if

  UT_Sun_in_south = ( RA - GMST0 - longitude ) / 15.0
  if (UT_Sun_in_south < 0) then
    UT_Sun_in_south=UT_Sun_in_south + 24
  end if

  LHA= sind(h) - (sind(latitude) * sind(Dec))/(cosd(latitude) * cosd(Dec))
  if (LHA > -1 AND LHA < 1) then
    LHA=acosd(LHA)/15
  else
    SunTimes = "No sunrise,No sunset"
    exit function
  end if
  hour_rise=UT_Sun_in_south - LHA
  hour_set=UT_Sun_in_south + LHA
  min_rise=int((hour_rise-int(hour_rise)) * 60)
  min_set=int((hour_set-int(hour_set)) * 60)

  hour_rise=(int(hour_rise) + (TZ + isdst))
  hour_set=(int(hour_set) + (TZ + isdst))
  if (min_rise < 10) then
    min_rise = right("0000" & min_rise, 2)
  end if
  if (min_set < 10) then
    min_set = right("0000" & min_set, 2)
  end if

  SunTimes = hour_rise & ":" & min_rise & "," & hour_set & ":" & min_set
End Function

' Support Functions
Function sind(qqq)
sind = sin((qqq) * DEGRAD)
End Function

Function cosd(qqq)
cosd = cos((qqq) * DEGRAD)
End Function

Function tand(qqq)
tand = tan((qqq) * DEGRAD)
End Function

Function atand(qqq)
atand = (RADEG * atan(qqq))
End Function

Function asind(qqq)
asind = (RADEG * asin(qqq))
End Function

Function acosd(qqq)
acosd = (RADEG * acos(qqq))
End Function

Function atan2d (qqq, qqq1)
atan2d = (RADEG * atan2(qqq, qqq1))
End Function

Function rev(qqq)
Dim x
x = (qqq - int(qqq/360.0) * 360.0)
if (x <= 0) then
x = x + 360
end if
rev = x
End Function

Function atan2(ys,xs)
' from the Draw Arrows tip
' @ http://www.devx.com/upload/free/features/VBPJ/techtips/techtips11.pdf
' by —Jim Deutch, Syracuse, New York

Dim theta
If xs <> 0 Then
  theta = Atn(ys / xs)
  If xs < 0 Then
   theta = theta + pi
  End If
Else
  If ys < 0 Then
   theta = 3 * pi / 2 '90
  Else
   theta = pi / 2 '270
  End If
End If
atan2 = theta
End Function

function acos(x)
acos = Atn(-X / Sqrt(-X * X + 1)) + 2 * Atn(1)
end function

function sqrt(x)
if x > 0 then
sqrt = Sqr(x)
else
sqrt = 0
end if
end function