VBA Code to calculate Sunrise and Sunset times in Auckland, NZ. Will need a bit of a tweak to Parameterize and Pythonorize it but it is all proven in the field.
[code]Sub Main
Dim vSunRiseTime As String
Dim vSunSetTime As String
Dim Season As String
Dim vTime As Date
vTime = Now
If GetTimeZone = "New Zealand Daylight Time" Then
Season = "Summer"
Else
Season = "Winter"
End If
vSunRiseTime = SunRiseSet("SunRise", 4) '4 minutes past sun rise
vSunSetTime = SunRiseSet("SunSet", 30) '30 minutes past sunset
End Sub
Private Function SunRiseSet(vSunRiseOrSunSet As String, vMinuteOffset As Integer) As String
Dim vTimeOffset As Double
Dim vLatitude As Double
Dim vLongitude As Double
Dim vZenith As Double
Dim vTimeZone As Integer 'Number of hours from UTC
Dim pi As Double
Dim vYear As Integer
Dim vDays As Integer
Dim lngHour As Double 'Longitudinal hour
Dim SinDec, CosDec As Double
Dim T As Double 'Local mean time of rising / setting
Dim Mdeg As Double
Dim Ldeg As Double
Dim RA As Double
Dim LQuadrant, RAQuadrant As Integer
Dim COSh As Double
Dim H As Double
Dim UT As Double
Dim LocalT As Double
Dim vHour As Integer
Dim vMinute As Integer
pi = 4 * Atn(1)
vLongitude = 174.766667 'Degrees: Auckland
vLatitude = -36.866667 'Degrees: Auckland
vZenith = 90.833 'Sun's zenith for sunrise/sunset
If Season = "Summer" Then
vTimeZone = 13
Else
vTimeZone = 12
End If
'Find the number of days so far this year
vYear = DatePart("yyyy", vTime)
vDays = DateDiff("d", DateSerial(vYear, 1, 1), vTime) + 1
'Convert longitude to hour value and calculate the approximate time
lngHour = vLongitude / (360 / 24)
'Approximate rising time in decimal days
If vSunRiseOrSunSet = "SunRise" Then T = vDays + ((6 - lngHour) / 24)
'Approximate setting time in decimal days
If vSunRiseOrSunSet = "SunSet" Then T = vDays + ((18 - lngHour) / 24)
'Calculate the Sun's mean anomaly
Mdeg = (0.9856 * T) - 3.289 'Degrees
Ldeg = Limit360(Mdeg + (1.916 * Sin(Radians(Mdeg))) + (0.02 * Sin(2 * Radians(Mdeg))) + 282.634) 'Degrees
RA = Limit360((180 / pi) * Atn(0.91764 * Tan(Radians(Ldeg))))
LQuadrant = (Floor(Ldeg / 90)) * 90
RAQuadrant = (Floor(RA / 90)) * 90
RA = RA + (LQuadrant - RAQuadrant)
RA = RA / 15
SinDec = 0.39782 * Sin(Radians(Ldeg))
CosDec = Cos(ArcSin(SinDec))
COSh = (Cos(Radians(vZenith)) - (SinDec * Sin(Radians(vLatitude)))) / (CosDec * Cos(Radians(vLatitude)))
’ if if rising time is desired:
If vSunRiseOrSunSet = “SunRise” Then H = 360 - ((180 / pi) * ArcCos(COSh))
’ if setting time is desired:
If vSunRiseOrSunSet = “SunSet” Then H = (180 / pi) * ArcCos(COSh)
H = H / 15
T = H + RA - (0.06571 * T) - 6.622
UT = Limit24(T - lngHour)
LocalT = Limit24(UT + vTimeZone) 'Decimal hours
vHour = Floor(LocalT)
vMinute = (LocalT - vHour) * 60
'Apply offset
vMinute = vMinute + vMinuteOffset
If vMinute >= 60 Then
vMinute = vMinute - 60
vHour = vHour + 1
End If
SunRiseSet = Format(vHour, "00") & ":" & Format(vMinute, "00")
End Function
Private Function Limit24(vX As Variant) As Double
Do
If vX > 24 Then vX = vX - 24
Loop Until vX < 24
Do
If vX < 0 Then vX = vX + 24
Loop Until vX >= 0
Limit24 = vX
End Function
Private Function Limit360(vX As Variant) As Double
Do
If vX > 360 Then vX = vX - 360
Loop Until vX < 360
Do
If vX < 0 Then vX = vX + 360
Loop Until vX >= 0
Limit360 = vX
End Function
Private Function Radians(vX As Variant) As Double
Dim pi As Double
pi = 4 * Atn(1)
Radians = vX * pi / 180
End Function
Private Function ArcCos(vX As Variant) As Variant
'Radians in, radians out
'Arccos(X) = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)
ArcCos = Atn((-1 * vX) / Sqr(-1 * vX * vX + 1)) + (2 * Atn(1))
End Function
Private Function ArcSin(vX As Variant) As Variant
'Arcsin(X) = Atn(X / Sqr(-X * X + 1))
ArcSin = Atn(vX / Sqr(-1 * vX * vX + 1))
End Function
Public Function Ceiling(ByVal X As Double, Optional ByVal Factor As Double = 1) As Double
’ X is the value you want to round
’ is the multiple to which you want to round
Ceiling = (Int(X / Factor) - (X / Factor - Int(X / Factor) > 0)) * Factor
End Function
Public Function Floor(ByVal X As Double, Optional ByVal Factor As Double = 1) As Double
’ X is the value you want to round
’ is the multiple to which you want to round
Floor = Int(X / Factor) * Factor
End Function
Private Type SYSTEMTIME
wYear As Integer
wMonth As Integer
wDayOfWeek As Integer
wDay As Integer
wHour As Integer
wMinute As Integer
wSecond As Integer
wMilliseconds As Integer
End Type
Private Type TIME_ZONE_INFORMATION
Bias As Long
StandardName(0 To 31) As Integer
StandardDate As SYSTEMTIME
StandardBias As Long
DaylightName(0 To 31) As Integer
DaylightDate As SYSTEMTIME
DaylightBias As Long
End Type
Private Declare Function GetTimeZoneInformation Lib “kernel32.dll” (lpTimeZoneInformation As TIME_ZONE_INFORMATION) As Long
Public Function GetTimeZone(Optional ByRef strTZName As String) As String
Dim objTimeZone As TIME_ZONE_INFORMATION
Dim lngResult As Long
Dim i As Long
lngResult = GetTimeZoneInformation&(objTimeZone)
Select Case lngResult
Case 0&, 1& 'use standard time
GetTimeZone = -(objTimeZone.Bias + objTimeZone.StandardBias) 'into minutes
For i = 0 To 31
If objTimeZone.StandardName(i) = 0 Then Exit For
strTZName = strTZName & Chr(objTimeZone.StandardName(i))
Next
Case 2& 'use daylight savings time
GetTimeZone = -(objTimeZone.Bias + objTimeZone.DaylightBias) 'into minutes
For i = 0 To 31
If objTimeZone.DaylightName(i) = 0 Then Exit For
strTZName = strTZName & Chr(objTimeZone.DaylightName(i))
Next
End Select
GetTimeZone = strTZName
End Function[/code]
