Here is an implementation of the Sunrise/Sunset Algorithm in VBA.
Official times of sunrise og sunset are calculated with variables depending on location given in latitude and longitude and corrected for difference between UCT/GMT and local time. Tested in Excel 2007.
Source of pseudocode algorithm: Almanac for Computers, 1990; Nautical Almanac Office; United States Naval Observatory; Washington, DC 20392
Inputs: day (integer), month (integer), year (double), latitude (double), longitude (double), localOffset (difference between local time and UCT/GMT), status (return time of official sun rise=0; return time of official sun set=1)
day, month and year are return values of standard Excel functions DAY, MONTH and YEAR applied to a cell formatted as an date.
My own location: 57.07 north and 9.94 east (eastern coordinates are positive, western negative) This year the sun will rise at 5:49:04 AM and set at 9:00:03 PM at my home address on my birthday.
Option Explicit
Public Const cPI As Double = 3.14159265358979
Public Function IntervalCorrection(a, b) Dim Temp As Double Temp = a If Temp < 0 Then Temp = Temp + b End If If Temp > b Then Temp = Temp – b End If IntervalCorrection = Temp End Function
‘Implementation of trigonometric functions taking degrees as argument instead of radians Public Function cosDeg(degree As Double) cosDeg = Cos(degree * cPI / 180#) End Function Public Function sinDeg(degree As Double) sinDeg = Sin(degree * cPI / 180#) End Function Public Function tanDeg(degree As Double) tanDeg = Tan(degree * cPI / 180#) End Function ‘Implementation of trigonometric functions returning degrees instead of radians Public Function acosDeg(arg As Double) acosDeg = 180# * Application.WorksheetFunction.Acos(arg) / cPI End Function Public Function asinDeg(arg As Double) asinDeg = 180# * Application.WorksheetFunction.Asin(arg) / cPI End Function Public Function atanDeg(arg As Double) atanDeg = 180# * Atn(arg) / cPI End Function
Public Function suntimes(day As Integer, month As Integer, year As Integer, latitude As Double, longitude As Double, localOffset As Integer, state As Double) ‘The day of the year and associated intermediate variables Dim N As Double Dim N1 As Double Dim N2 As Double Dim N2A As Double Dim N3 As Double ‘Hour value and approximate time Dim lngHour As Double Dim t As Double ‘Suns mean anomaly Dim M As Double ‘Suns true longitude and intermediate variables Dim L As Double Dim L1 As Double Dim L2 As Double ‘Suns right ascension Dim RA As Double ‘Intermediate variables in calculation of ascension value Dim Lquadrant As Double Dim RAquadrant As Double ‘Suns declination Dim sinDec As Double Dim cosDec As Double ‘Suns local hour value and intermediate values Dim H As Double Dim cosH As Double Dim cosH1 As Double Dim cosH2 As Double Dim cosH3 As Double ‘Local mean time of rising Dim LT As Double ‘Coordinated Universal Time of rising Dim UT As Double ‘Definition of zenith – official 90.0 degrees, civil 96 50’ degrees, nautical 102 degrees, astronomical 108 degrees Dim zenith As Double zenith = 90#
‘Calculate the day of the year N1 = WorksheetFunction.Floor(275# * month / 9#, 1) N2 = WorksheetFunction.Floor((month + 9#) / 12#, 1) N2A = WorksheetFunction.Floor(year / 4#, 1) N3 = 1 + WorksheetFunction.Floor((year – 4# * N2A + 2#) / 3#, 1) N = N1 – (N2 * N3) + day – 30# Debug.Print “The value of variable N is: ” & N ‘Convert the longitude to hour value and calculate an approximate time lngHour = longitude / 15# t = N + ((6# – lngHour) / 24#) If (state > 0) Then t = N + ((18# – lngHour) / 24#) End If
‘Calculate the Sun’s mean anomaly M = (0.9856 * t) – 3.289 Debug.Print “The value of variable M is: ” & M ‘Calculate the Suns true longitude L1 = 1.916 * sinDeg(M) L2 = 0.02 * sinDeg(2# * M) L = M + L1 + L2 + 282.634 L = IntervalCorrection(L, 360)
‘Calculate the Suns right ascension RA = atanDeg(0.91764 * tanDeg(L)) RA = IntervalCorrection(RA, 360)
‘Right ascension needs to be in the same quadrant as L Lquadrant = (Application.WorksheetFunction.Floor(L / 90#, 1)) * 90# RAquadrant = (Application.WorksheetFunction.Floor(RA / 90#, 1)) * 90# RA = RA + (Lquadrant – RAquadrant)
‘Right ascension value needs to be converted into hours RA = RA / 15#
‘Calculate the Suns declination sinDec = 0.39782 * sinDeg(L) cosDec = cosDeg(asinDeg(sinDec))
‘Calculate the Suns local hour angle cosH1 = cosDeg(zenith) cosH2 = sinDec * sinDeg(latitude) cosH3 = cosDec * cosDeg(latitude) cosH = (cosH1 – cosH2) / cosH3
If (cosH > 1) Then MsgBox (“The sun do not rise”) Exit Function End If
If (cosH < -1) Then MsgBox (“The sun do not set”) Exit Function End If
‘Finish calculating H and convert into hours H = (360 – acosDeg(cosH)) / 15# If (state > 0) Then H = acosDeg(cosH) / 15# End If
‘Calculate local mean time of rising LT = H + RA – (0.06571 * t) – 6.622 ‘Adjust back to UTC UT = LT – lngHour UT = IntervalCorrection(UT, 24) Debug.Print “The value of variable cPI is: ” & cPI Debug.Print “The value of variable UT is: ” & UT Debug.Print “The value of variable localOffset is: ” & localOffset ‘Return value suntimes = UT + localOffset End Function