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

Comments

Popular Posts