Sub Class_Globals
Private PI As Double = 3.141592653589793
Private EarthRadius As Double = 6371.0 ' Earth's radius in kilometers
Private AtmosphericRefraction As Double = 34.0 / 60.0 ' in degrees, typically 34 arcminute
Private fixMinutes As Double
End Sub
Public Sub Initialize
End Sub
Sub CalculateSunriseSunset(lat As Double, lon As Double, elev As Double, tz As Double, year As Int, month As Int, day As Int, isSunrise As Boolean) As Double
' Convert the date to Julian date
Dim Jday As Double = JulianDate(year, month, day)
Dim Jcentury As Double = (Jday - 2451545.0) / 36525.0
' Calculate the sun's mean anomaly
Dim M As Double = 357.52911 + Jcentury * (35999.05029 - 0.0001537 * Jcentury)
M = Bit.FMod(M,360)
' Calculate the sun's true longitude
Dim C As Double = (1.914602 - Jcentury * (0.004817 + 0.000014 * Jcentury)) * Sin(DegToRad(M)) + (0.019993 - 0.000101 * Jcentury) * Sin(DegToRad(2 * M)) + 0.000289 * Sin(DegToRad(3 * M))
Dim Lsun As Double = Bit.FMod((280.46646 + Jcentury * (36000.76983 + 0.0003032 * Jcentury) + C),360)
' Calculate the sun's right ascension
Dim RA As Double = RadToDeg(ATan(0.91764 * Tan(DegToRad(Lsun))))
RA = Bit.FMod(RA + 360,360)
' Right ascension value needs to be in the same quadrant as Lsun
Dim Lquadrant As Int = Floor(Lsun / 90) * 90
Dim RAquadrant As Int = Floor(RA / 90) * 90
RA = RA + (Lquadrant - RAquadrant)
RA = RA / 15 ' Convert to hours
' Calculate the sun's declination
Dim sinDec As Double = 0.39782 * Sin(DegToRad(Lsun))
Dim cosDec As Double = Cos(ASin(sinDec))
' Calculate the sun's local hour angle
Dim zenith As Double = 90.833 + AtmosphericRefraction + AdjustForElevation(elev)
Dim cosH As Double = (Cos(DegToRad(zenith)) - sinDec * Sin(DegToRad(lat))) / (cosDec * Cos(DegToRad(lat)))
' Ensure cosH is within the range of -1 to 1
If cosH > 1 Then
Return -1 ' The sun never rises on this location on the specified date
Else If cosH < -1 Then
Return -1 ' The sun never sets on this location on the specified date
End If
Dim H As Double
If isSunrise Then
H = 360 - RadToDeg(ACos(cosH))
Else
H = RadToDeg(ACos(cosH))
End If
H = H / 15 ' Convert to hours
' Calculate local mean time of rising/setting
Dim T As Double = H + RA - (0.06571 * (Jday - 2451545.0)) - 6.622
' Adjust back to UTC
Dim UT As Double = T - (lon / 15)
UT = Bit.FMod((UT + 24), 24) ' Ensure UT is within 0-24 range
' Convert UTC to local time
Dim localT As Double = UT + tz
localT = Bit.FMod(localT + 24,24) ' Ensure localT is within 0-24 range
If isSunrise Then 'fix for sunset
fixMinutes = 1.25/60 '
localT = localT - fixMinutes
Else
fixMinutes = 7/60 '
localT = localT - fixMinutes
End If
If localT < 0 Then localT = localT + 24 'fix if is negative value
Return localT
End Sub
public Sub JulianDate(year As Int, month As Int, day As Int) As Double
If month <= 2 Then
year = year - 1
month = month + 12
End If
Dim A As Int = Floor(year / 100)
Dim B As Int = 2 - A + Floor(A / 4)
Return Floor(365.25 * (year + 4716)) + Floor(30.6001 * (month + 1)) + day + B - 1524.5
End Sub
Sub DegToRad(deg As Double) As Double
Return deg * PI / 180
End Sub
Sub RadToDeg(rad As Double) As Double
Return rad * 180 / PI
End Sub
public Sub ConvertToTime(d As Double) As String
Dim hour As Int = Floor(d)
Dim minute As Int = Floor((d - hour) * 60)
Dim seconds As Int = (((d - hour) * 60) - Floor((d - hour) * 60)) * 60
Return NumberFormat2(hour, 2, 0, 0, False) & ":" & NumberFormat2(minute, 2, 0, 0, False) & ":" & NumberFormat2(seconds, 2, 0, 0, False)
End Sub
public Sub TimeToDouble(TimeString As String) As Double
Dim parts() As String = Regex.Split(":", TimeString)
Dim hours As Int = parts(0)
Dim minutes As Int = parts(1)
Dim seconds As Int = 0
If parts.Length = 3 Then seconds = parts(2)
Dim totalSeconds As Int = (hours * 3600) + (minutes * 60) + seconds
Dim totalSecondsInDay As Int = 24 * 3600
Return totalSeconds / totalSecondsInDay
End Sub
public Sub AdjustForElevation(elev As Double) As Double
Dim angle As Double = RadToDeg(ACos(EarthRadius / (EarthRadius + elev / 1000)))
Return angle
End Sub
