'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