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]