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

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