001  <%
002   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
003   '::: :::
004   '::: These functions calculate sunrise and sunset times for any given :::
005   '::: latitude and longitude. They may also be used to calculate such :::
006   '::: things as astronomical twilight, nautical twilight and civil :::
007   '::: twilight. :::
008   '::: :::
009   '::: SPECIAL NOTES: This code is valid for dates from 1901 to 2099, and :::
010   '::: will not calculate sunrise/set times for latitudes :::
011   '::: above/below 63/-63 degrees. :::
012   '::: :::
013   '::: This code is based on the work of several others, including Jean :::
014   '::: Meeus, Todd Guillory, Christophe David, Kieth Burnett and Roger W. :::
015   '::: Sinnott (credit where due!) :::
016   '::: :::
017   '::: Converted to VBScript and cleaned-up by Mike Shaffer. :::
018   '::: :::
019   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
020  
021   const pi = 3.14159265358979323846
022  
023   degrees = 180 / pi
024   radians = pi / 180
025  
026   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
027   '::: Returns an angle in range of 0 to (2 * pi) :::
028   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
029   function GetRange (x)
030  
031   dim temp1
032   dim temp2
033  
034   temp1 = x / (2 * pi)
035   temp2 = (2 * pi) * (temp1 - fix(temp1))
036   IF temp2 < 0 THEN
037   temp2 = (2 * pi) + temp2
038   End If
039  
040   GetRange = temp2
041  
042   END function
043  
044  
045   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
046   '::: Returns 24 hour time from decimal time :::
047   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
048   Function GetMilitaryTime(DecimalTime, GMTOffset)
049  
050   dim temp1
051   dim temp2
052  
053   ' Handle 24-hour time wrap
054   if decimaltime+gmtoffset < 0 then decimaltime = decimaltime + 24
055   if decimaltime+gmtoffset > 24 then decimaltime = decimaltime - 24
056  
057   temp1 = ABS(DecimalTime + GMTOffset)
058   temp2 = INT(temp1)
059   temp1 = 60 * (temp1 - temp2)
060   temp1 = right("0000" & CSTR(INT(temp2 * 100 + temp1 + .5)), 4)
061  
062   GetMilitaryTime = left(temp1, 2) & ":" & right(temp1, 2)
063  
064   END Function
065  
066  
067   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
068   '::: This routine does all the real work :::
069   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
070   function GetSunRiseSet(latitude, ByVal longitude, ZoneRelativeGMT, RiseOrSet, Year, Month, Day)
071   if abs(latitude) > 63 then
072   GetSunRiseSet = "{invalid latitude}"
073   Exit Function
074   end if
075  
076   ' An altitude of -0.833 is generally accepted as the angle of
077   ' the sun at which sunrise/sunset occurs. It is not exactly
078   ' zero because of refraction effects of the earth's atmosphere.
079   Altitude = -0.833
080  
081   select case ucase(RiseOrSet)
082   case "S"
083   RS = -1
084   case else
085   RS = 1
086   end select
087  
088   Ephem2000Day = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5
089  
090   utold = pi
091   utnew = 0
092   sinalt = cdbl(SIN(altitude * radians)) ' solar altitude
093   sinphi = cdbl(SIN(latitude * radians)) ' viewer's latitude
094   cosphi = cdbl(COS(latitude * radians)) '
095   longitude = cdbl(longitude * radians) ' viewer's longitude
096  
097   Err.Clear
098   On Error Resume Next
099  
100   DO WHILE (ABS(utold - utnew) > .001) and (ct < 35)
101  
102   ct = ct + 1
103  
104   utold = utnew
105   days = Ephem2000Day + utold / (2 * pi)
106   t = days / 36525
107   ' These 'magic' numbers are orbital elements of the sun, and should not be changed
108   L = GetRange(4.8949504201433 + 628.331969753199 * t)
109   G = GetRange(6.2400408 + 628.3019501 * t)
110   ec = .033423 * SIN(G) + .00034907 * SIN(2. * G)
111   lambda = L + ec
112   E = -1 * ec + .0430398 * SIN(2. * lambda) - .00092502 * SIN(4. * lambda)
113   obl = .409093 - .0002269 * t
114  
115   ' Obtain ASIN of (SIN(obl) * SIN(lambda))
116   delta = SIN(obl) * SIN(lambda)
117   delta = ATN(delta / (sqr(1 - delta * delta)))
118  
119   GHA = utold - pi + E
120   cosc = (sinalt - sinphi * SIN(delta)) / (cosphi * COS(delta))
121   SELECT CASE cosc
122   CASE cosc > 1
123   correction = 0
124   CASE cosc < -1
125   correction = pi
126   CASE ELSE
127   correction = atn((sqr(1 - cosc * cosc)) / cosc)
128   END SELECT
129  
130   utnew = GetRange(utold - (GHA + longitude + RS * correction))
131  
132   LOOP
133  
134   If err = 0 then
135   GetSunRiseSet = GetMilitaryTime(utnew * degrees / 15, ZoneRelativeGMT)
136   else
137   GetSunRiseSet = "{err}"
138   end if
139  
140   end function
141  
142  
143  
144  
145  
146   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
147   '::: Here is test harness. This code should be removed for production. :::
148   '::: :::
149   '::: This example code shows today's sunrise time for Dallas :::
150   ':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
151  
152   y = year(now)
153   m = month(now)
154   d = day(now)
155  
156   ' Set these to the latitude/longitude of the location
157   MyLatitude = 32.9697
158   MyLongitude = -96.80322
159  
160   ' Set this to your offset from GMT (e.g. for Dallas is -6)
161   ' NOTE: The routine does NOT handle switches to/from daylight savings
162   ' time, so beware!
163   MyTimeZone = -6
164  
165   ' Note:Set RiseOrSet to "R" for sunrise, "S" for sunset
166   RiseOrSet = "R"
167   ttt = GetSunRiseSet(MyLatitude, MyLongitude, MyTimeZone, RiseOrSet, y, m, d)
168  
169   RiseOrSet = "S"
170   sss = GetSunRiseSet(MyLatitude, MyLongitude, MyTimeZone, RiseOrSet, y, m, d)
171  
172   response.write "<html><body>Sunrise for Dallas on " & _
173   datevalue(now) & " occurs at: " & ttt & _
174   " (24-hour time)<P>Sunset occurs at " & _
175   sss & "</body></html>"
176  
177  %>