PaulSadowski.com 

Windows Scripting Host
 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 separated string indicating the sunrise and sunset times (24 hour clock), respectively. 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 strOu
t

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

 

 

¬© 2003 by Paul R. Sadowski   
All Rights Reserved. Used By Permission.  
Comments to: Scripting