B4J Code Snippet sunrise and sunset sub

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

public Sub AdjustForElevation(elev As Double) As Double
    Dim angle As Double = RadToDeg(ACos(EarthRadius / (EarthRadius + elev / 1000)))
    Return angle
End Sub

dfo = rads*(0.5*SunDia + AirRefr)