Solar Position UDF - Debug

SirAlonne · April 10, 2020

Hi,

I recently found this UDF and it gives pretty accurate sunrise and sunset times compared to google (no difference), but there are some instances when a repetition occurs..

2020. 03. 30-31 produces the same sunrise and sunset time (where normally should be a difference between the two days)
("_atan2" function was missing from the original)

Maybe someone could help me figure out what causes this.

#include-once
#include <Array.au3>
#include <Math.au3>
#include <Date.au3>

Global Const $Pi = 4 * ATan(1)
Global Const $MonthLength[12] = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
Global Const $MonthFullName[12] = ["January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"]
Global Const $MonthAbrev[12] = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]


; #FUNCTION# ====================================================================================================================
; Name...........: _GetSolarAstronimicalData
; Description ...: Calculate the position of the sun relative to a position on earth for a given time/date.
; Syntax.........: _GetSolarAstronimicalData ( $Latitude, $Longitude, $Month, $MDay, $Year, $TimeZone, $Hour, $Minute, $Second, $DST )
; Parameters ....: $Latitude     - The Latitude portion of the earth position to calculate the sun's position relative to (float)
;                  $Longitude     - The Latitude portion of the earth position to calculate the sun's position relative to (float)
;                  $Month     - The month of the year you are calculating for (numeric-string: "1-12")
;                  $MDay     - The day of the month you are calculating for (numeric-string: "1-31")
;                  $Year     - The year you are calculating for (numeric-string: "-2000 to 3000")
;                  $TimeZone     - The time zone you are calculating for (numeric-string: "-11 to 12")
;                  $Hour     - Optional!  The hour you are calculating for in 24-hrs (numeric-string: "0-23")
;                  $Minute     - Optional!  The minute you are calculating for (numeric-string: "00-59")
;                  $Second    - Optional!  The second you are calculating for (numeric-string: "00-59")
;                  $DST     - Optional!  Is Daylight Saving's Time in effect? (boolean)
; Return values .: Success - Returns an array with the following values:
;                            0 = Element count. If count = 3, then only items 1-3 available.  If count = 6, then all elements available)
;                            1 = Sunrise time "07:12" 24-hr time, or "07:12 Jul 17" if position is near international date line
;                            2 = Solar Noon "13:16:46" 24-hr time with seconds.  Higest point of sun in sky.
;                            3 = Sunset time "19:22" 24-hr time, or "19:22 Jul 19" if position is near international date line
;                            4 = Sun Azimuth in degrees relative to computed position (0deg = true north, 180deg = true south, 90deg = East, 270 = West)
;                            5 = Sun Elevation is degrees relative to computed position at sea-level
;                            6 = Ambient light disposition (possibilities: Day, Civil Twilight, Nautical Twilight, Astronomical Twilight, Night)
; Author ........: Tim Strommen (tim292stro)
; Modified.......:
; Remarks .......: Again, special thanks to NOAA for allowing re-use of their code!!  See: "http://www.esrl.noaa.gov/"
; Related .......:
; Link ..........:
; Example .......: Yes
; ===============================================================================================================================








;GoogleMaps | UK Greenwich
Global $Latitude = 51.4874009
Global $Longitude = -0.012965
Global $TimeZone = 0
Global $Test, $DST = 0

For $x = 30 To 31 Step 1
    $Test = _GetSolarAstronimicalData($Latitude, $Longitude, 3, $x, 2020, $TimeZone)
    _ArrayDisplay($Test, "(Leap Year) 2020 March " & $x)
Next

















Func _GetSolarAstronimicalData($AstroLat = "", $AstroLong = "", $AstroMonth = "", $AstroDay = "", $AstroYear = "", $AstroTimeZone = "", $AstroHour = "", $AstroMinute = "", $AstroSecond = "", $AstroDST = False)
    Local $AstroBeginAT, $AstroBeginNT, $AstroBeginCT, $AstroSunrise, $AstroSolarNoon, $AstroSunset, $AstroEndCT, $AstroEndNT, $AstroEndAT, $AstroSolarElevation, $AstroSolarAzimuth
    If ($AstroLat = "") Or ($AstroLong = "") Then
        ; Return current Greenwich, UK basic times with current solar position
        $AstroLat = "51.48"
        $AstroLong = "0.000000001"
        $AstroTimeZone = "0"
        $AstroMonth = @MON
        $AstroDay = @MDAY
        $AstroYear = @YEAR
        $AstroHour = @HOUR
        $AstroMinute = @MIN
        $AstroSecond = @SEC
        $RetunedArray = _SolarCalculate($AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST, True)
        Return $RetunedArray
    ElseIf $AstroHour = "" Then
        ; Just return the specified day's basic times, no current solar position
        $AstroHour = "12"
        $AstroMinute = "00"
        $AstroSecond = "00"
        $RetunedArray = _SolarCalculate($AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST, False)
        Return $RetunedArray
    Else
        ; Return both the basic times, plus the current solar position
        $RetunedArray = _SolarCalculate($AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST, True)
        Return $RetunedArray
    EndIf
EndFunc   ;==>_GetSolarAstronimicalData


Func _CalcTimeJulianCent($jd = 0)
    Local $Time = ($jd - 2451545.0) / 36525.0
    Return $Time
EndFunc   ;==>_CalcTimeJulianCent


Func _CalcJDFromJulianCent($t = 0)
    Local $jd = $t * 36525.0 + 2451545.0
    Return $jd
EndFunc   ;==>_CalcJDFromJulianCent


Func _isLeapYear($yr = 0)
    Return ((Mod($yr, 4) And Mod($yr, 100)) Or (Mod($yr, 400) = 0))
EndFunc   ;==>_isLeapYear


Func _calcDoyFromJD($jd = 0)
    Local $z = Floor($jd + 0.5)
    Local $f = ($jd + 0.5) - $z
    Local $A, $alpha, $B, $C, $D, $doy, $E, $day, $month, $year, $K
    If $z < 2299161 Then
        $A = $z
    Else
        $alpha = Floor(($z - 1867216.25) / 36524.25)
        $A = $z + 1 + $alpha - Floor($alpha / 4)
    EndIf
    $B = $A + 1524
    $C = Floor(($B - 122.1) / 365.25)
    $D = Floor(365.25 * $C)
    $E = Floor(($B - $D) / 30.6001)
    $day = $B - $D - Floor(30.6001 * $E) + $f
    If $E < 14 Then
        $month = $E - 1
    Else
        $month = $E - 13
    EndIf
    If $month > 2 Then
        $year = $C - 4716
    Else
        $year = $C - 4715
    EndIf
    If _isLeapYear($year) Then
        $K = 1
    Else
        $K = 2
    EndIf
    $doy = Floor((275 * $month) / 9) - $K * Floor(($month + 9) / 12) + $day - 30
    Return $doy
EndFunc   ;==>_calcDoyFromJD


Func _radToDeg($AngleRad = 0)
    Return 180.0 * $AngleRad / $Pi
EndFunc   ;==>_radToDeg


Func _degToRad($AngleDeg = 0)
    Return $Pi * $AngleDeg / 180.0
EndFunc   ;==>_degToRad


Func _atan2($ys, $xs)
    Local $theta
    If $xs <> 0 Then
        $theta = ATan($ys / $xs)
        If $xs < 0 Then
            $theta = $theta + $Pi
        EndIf
    Else
        If $ys < 0 Then
            $theta = 3 * $Pi / 2 ;90
        Else
            $theta = $Pi / 2 ;270
        EndIf
    EndIf
    $atan2 = $theta
    Return $atan2
EndFunc   ;==>_atan2


Func _calcGeomMeanLongSun($t = 0)
    Local $L0 = 280.46646 + $t * (36000.76983 + $t * (0.0003032))
    While $L0 > 360.0
        $L0 -= 360.0
    WEnd
    While $L0 < 0.0
        $L0 += 360.0
    WEnd
    Return $L0 ; in degrees
EndFunc   ;==>_calcGeomMeanLongSun


Func _calcGeomMeanAnomalySun($t = 0)
    Local $M = 357.52911 + $t * (35999.05029 - 0.0001537 * $t)
    Return $M  ; in degrees
EndFunc   ;==>_calcGeomMeanAnomalySun


Func _calcEccentricityEarthOrbit($t = 0)
    Local $E = 0.016708634 - $t * (0.000042037 + 0.0000001267 * $t)
    Return $E  ; unitless
EndFunc   ;==>_calcEccentricityEarthOrbit


Func _calcSunEqOfCenter($t = 0)
    Local $M = _calcGeomMeanAnomalySun($t)
    Local $mrad = _degToRad($M)
    Local $sinm = Sin($mrad)
    Local $sin2m = Sin($mrad + $mrad)
    Local $sin3m = Sin($mrad + $mrad + $mrad)
    Local $C = $sinm * (1.914602 - $t * (0.004817 + 0.000014 * $t)) + $sin2m * (0.019993 - 0.000101 * $t) + $sin3m * 0.000289
    Return $C
EndFunc   ;==>_calcSunEqOfCenter


Func _calcSunTrueLong($t = 0)
    Local $lo = _calcGeomMeanLongSun($t)
    Local $C = _calcSunEqOfCenter($t)
    Local $O = $lo + $C
    Return $O  ; in degrees
EndFunc   ;==>_calcSunTrueLong


Func _calcSunTrueAnomaly($t = 0)
    Local $M = _calcGeomMeanAnomalySun($t)
    Local $C = _calcSunEqOfCenter($t)
    Local $v = $M + $C
    Return $v
EndFunc   ;==>_calcSunTrueAnomaly


Func _calcSunRadVector($t = 0)
    Local $v = _calcSunTrueAnomaly($t)
    Local $E = _calcEccentricityEarthOrbit($t)
    Local $R = (1.000001018 * (1 - $E * $E)) / (1 + $E * Cos(_degToRad($v)))
    Return $R  ; in AUs
EndFunc   ;==>_calcSunRadVector


Func _calcSunApparentLong($t = 0)
    Local $O = _calcSunTrueLong($t)
    Local $omega = 125.04 - 1934.136 * $t
    Local $lambda = $O - 0.00569 - 0.00478 * Sin(_degToRad($omega))
    Return $lambda  ; in degrees
EndFunc   ;==>_calcSunApparentLong


Func _calcMeanObliquityOfEcliptic($t = 0)
    Local $seconds = 21.448 - $t * (46.8150 + $t * (0.00059 - $t * (0.001813)))
    Local $e0 = 23.0 + (26.0 + ($seconds / 60.0)) / 60.0
    Return $e0  ; in degrees
EndFunc   ;==>_calcMeanObliquityOfEcliptic


Func _calcObliquityCorrection($t = 0)
    Local $e0 = _calcMeanObliquityOfEcliptic($t)
    Local $omega = 125.04 - 1934.136 * $t
    Local $E = $e0 + 0.00256 * Cos(_degToRad($omega))
    Return $E
EndFunc   ;==>_calcObliquityCorrection


Func _calcSunRtAscension($t = 0)
    Local $E = _calcObliquityCorrection($t)
    Local $lambda = _calcSunApparentLong($t)
    Local $tananum = (Cos(_degToRad($E)) * Sin(_degToRad($lambda)))
    Local $tanadenom = (Cos(_degToRad($lambda)))
    Local $alpha = _radToDeg(_atan2($tananum, $tanadenom))
    Return $alpha  ; in degrees
EndFunc   ;==>_calcSunRtAscension


Func _calcSunDeclination($t = 0)
    Local $E = _calcObliquityCorrection($t)
    Local $lambda = _calcSunApparentLong($t)
    Local $sint = Sin(_degToRad($E)) * Sin(_degToRad($lambda))
    Local $theta = _radToDeg(ASin($sint))
    Return $theta  ; in degrees
EndFunc   ;==>_calcSunDeclination


Func _calcEquationOfTime($t = 0)
    Local $epsilon = _calcObliquityCorrection($t)
    Local $L0 = _calcGeomMeanLongSun($t)
    Local $E = _calcEccentricityEarthOrbit($t)
    Local $M = _calcGeomMeanAnomalySun($t)
    Local $y = Tan(_degToRad($epsilon) / 2.0)
    $y *= $y
    Local $sin2l0 = Sin(2.0 * _degToRad($L0))
    Local $sinm = Sin(_degToRad($M))
    Local $cos2l0 = Cos(2.0 * _degToRad($L0))
    Local $sin4l0 = Sin(4.0 * _degToRad($L0))
    Local $sin2m = Sin(2.0 * _degToRad($M))
    Local $Etime = $y * $sin2l0 - 2.0 * $E * $sinm + 4.0 * $E * $y * $sinm * $cos2l0 - 0.5 * $y * $y * $sin4l0 - 1.25 * $E * $E * $sin2m
    Return _radToDeg($Etime) * 4.0     ; in minutes of time
EndFunc   ;==>_calcEquationOfTime

;SunRise = Horizon+0.8333
Func _calcHourAngleSunrise($lat, $solarDec)
    Local $latRad = _degToRad($lat)
    Local $sdRad = _degToRad($solarDec)
    Local $HAarg = (Cos(_degToRad(90.833)) / (Cos($latRad) * Cos($sdRad)) - Tan($latRad) * Tan($sdRad))
    Local $HA = ACos($HAarg)
    Return $HA  ;  in radians for morning (for evening, use -HA)
EndFunc   ;==>_calcHourAngleSunrise

;CivilTwilight = (Horizon+6) < SunCenter < (Horizon+0.8333)
Func _calcHourAngleCivilTwilight($lat, $solarDec)
    Local $latRad = _degToRad($lat)
    Local $sdRad = _degToRad($solarDec)
    Local $HAarg = (Cos(_degToRad(96.0)) / (Cos($latRad) * Cos($sdRad)) - Tan($latRad) * Tan($sdRad))
    Local $HA = ACos($HAarg)
    Return $HA  ;  in radians for morning (for evening, use -HA)
EndFunc   ;==>_calcHourAngleCivilTwilight

;NauticalTwilight = (Horizon+12) < SunCenter < (Horizon+6)
Func _calcHourAngleNauticalTwilight($lat, $solarDec)
    Local $latRad = _degToRad($lat)
    Local $sdRad = _degToRad($solarDec)
    Local $HAarg = (Cos(_degToRad(102.0)) / (Cos($latRad) * Cos($sdRad)) - Tan($latRad) * Tan($sdRad))
    Local $HA = ACos($HAarg)
    Return $HA  ;  in radians for morning (for evening, use -HA)
EndFunc   ;==>_calcHourAngleNauticalTwilight

;AstronimicalTwilight = (Horizon+18) < SunCenter < (Horizon+12)
Func _calcHourAngleAstronimicalTwilight($lat, $solarDec)
    Local $latRad = _degToRad($lat)
    Local $sdRad = _degToRad($solarDec)
    Local $HAarg = (Cos(_degToRad(108.0)) / (Cos($latRad) * Cos($sdRad)) - Tan($latRad) * Tan($sdRad))
    Local $HA = ACos($HAarg)
    Return $HA  ;  in radians for morning ( for evening, use -HA )
EndFunc   ;==>_calcHourAngleAstronimicalTwilight


Func _isNumber($inputval)
    Return (IsFloat($inputval) Or IsNumber($inputval) Or IsInt($inputval))
EndFunc   ;==>_isNumber


Func _zeroPad($n, $digits)
    $n = String($n)
    While StringLen($n) < $digits
        $n = "0" & $n
    WEnd
    Return $n
EndFunc   ;==>_zeroPad

; Global variable defined at top of file...
Func _getJD($CompJDMon = "", $CompJDDay = "", $CompJDYear = "")
    If ((_isLeapYear($CompJDYear)) And ($CompJDMon = 2)) Then
        If $CompJDDay > 29 Then
            $CompJDDay = 29
        EndIf
    Else
        If $CompJDDay > $MonthLength[$CompJDMon] Then
            $CompJDDay = $MonthLength[$CompJDMon]
        EndIf
    EndIf
    If $CompJDMon <= 2 Then
        $CompJDYear -= 1
        $CompJDMon += 12
    EndIf
    Local $A = Floor($CompJDYear / 100)
    Local $B = 2 - $A + Floor($A / 4)
    Local $jd = Floor(365.25 * ($CompJDYear + 4716)) + Floor(30.6001 * ($CompJDMon + 1)) + $CompJDDay + $B - 1524.5
    Return $jd
EndFunc   ;==>_getJD


Func _getTimeLocal($CompHour = "", $CompMin = "", $CompSec = "", $CompDST = False)
    If ($CompDST) Then
        $CompHour -= 1
    EndIf
    Local $mins = $CompHour * 60 + $CompMin + $CompSec / 60.0
    Return $mins
EndFunc   ;==>_getTimeLocal


Func _calcAzEl($output, $t, $localtime, $Latitude, $Longitude, $zone)
    Local $SolarReturnAzEl[1] = [0]
    Local $SolarLightStatus, $SolarEquationOfTime, $SolarDeclination
    Local $eqTime = _calcEquationOfTime($t)
    Local $theta = _calcSunDeclination($t)
    If $output Then
        $SolarEquationOfTime = Floor($eqTime * 100 + 0.5) / 100.0
        $SolarDeclination = Floor($theta * 100 + 0.5) / 100.0
    EndIf
    Local $solarTimeFix = $eqTime + 4.0 * $Longitude - 60.0 * $zone
    Local $earthRadVec = _calcSunRadVector($t)
    Local $trueSolarTime = $localtime + $solarTimeFix
    While $trueSolarTime > 1440
        $trueSolarTime -= 1440
    WEnd
    Local $hourAngle = $trueSolarTime / 4.0 - 180.0 ;
    If $hourAngle < -180 Then
        $hourAngle += 360.0
    EndIf
    Local $haRad = _degToRad($hourAngle)
    Local $csz = Sin(_degToRad($Latitude)) * Sin(_degToRad($theta)) + Cos(_degToRad($Latitude)) * Cos(_degToRad($theta)) * Cos($haRad)
    If $csz > 1.0 Then
        $csz = 1.0
    ElseIf $csz < -1.0 Then
        $csz = -1.0
    EndIf
    Local $zenith = _radToDeg(ACos($csz))
    Local $azDenom = (Cos(_degToRad($Latitude)) * Sin(_degToRad($zenith)))
    If Abs($azDenom) > 0.001 Then
        $azRad = ((Sin(_degToRad($Latitude)) * Cos(_degToRad($zenith))) - Sin(_degToRad($theta))) / $azDenom
        If Abs($azRad) > 1.0 Then
            If $azRad < 0 Then
                $azRad = -1.0
            Else
                $azRad = 1.0
            EndIf
        EndIf
        Local $azimuth = 180.0 - _radToDeg(ACos($azRad))
        If $hourAngle > 0.0 Then
            $azimuth = -$azimuth
        EndIf
    Else
        If $Latitude > 0.0 Then
            $azimuth = 180.0
        Else
            $azimuth = 0.0
        EndIf
    EndIf
    If $azimuth < 0.0 Then
        $azimuth += 360.0
    EndIf
    Local $exoatmElevation = 90.0 - $zenith
    ; Atmospheric Refraction correction
    If $exoatmElevation > 85.0 Then
        Local $refractionCorrection = 0.0
    Else
        Local $te = Tan(_degToRad($exoatmElevation))
        If $exoatmElevation > 5.0 Then
            Local $refractionCorrection = 58.1 / $te - 0.07 / ($te * $te * $te) + 0.000086 / ($te * $te * $te * $te * $te)
        ElseIf $exoatmElevation > -0.575 Then
            Local $refractionCorrection = 1735.0 + $exoatmElevation * (-518.2 + $exoatmElevation * (103.4 + $exoatmElevation * (-12.79 + $exoatmElevation * 0.711)))
        Else
            Local $refractionCorrection = -20.774 / $te
        EndIf
        $refractionCorrection = $refractionCorrection / 3600.0
    EndIf

    Local $solarZen = $zenith - $refractionCorrection

    If $solarZen > 108.0 Then
        $SolarLightStatus = "Night"
    ElseIf ((108.0 > $solarZen) And ($solarZen >= 102.0)) Then
        $SolarLightStatus = "Astronomical Twilight"
    ElseIf ((102.0 > $solarZen) And ($solarZen >= 96.0)) Then
        $SolarLightStatus = "Nautical Twilight"
    ElseIf ((96.0 > $solarZen) And ($solarZen >= 90.8333)) Then
        $SolarLightStatus = "Civil Twilight"
    Else
        $SolarLightStatus = "Day"
    EndIf

    _ArrayAdd($SolarReturnAzEl, Floor((($azimuth * 100) + 0.5) / 100.0))
    _ArrayAdd($SolarReturnAzEl, Floor((90.0 - $solarZen) * 100 + 0.5) / 100.0)
    _ArrayAdd($SolarReturnAzEl, $SolarLightStatus)
    _ArrayDelete($SolarReturnAzEl, 0)

    Return ($SolarReturnAzEl)
EndFunc   ;==>_calcAzEl


Func _calcSolNoon($jd, $Longitude, $TimeZone, $DST)
    Local $tnoon = _CalcTimeJulianCent($jd - $Longitude / 360.0)
    Local $eqTime = _calcEquationOfTime($tnoon)
    Local $solNoonOffset = 720.0 - ($Longitude * 4) - $eqTime    ; in minutes
    Local $newt = _CalcTimeJulianCent($jd + $solNoonOffset / 1440.0)
    $eqTime = _calcEquationOfTime($newt)
    $solNoonLocal = 720 - ($Longitude * 4) - $eqTime + ($TimeZone * 60.0)      ; in minutes
    If $DST Then $solNoonLocal += 60.0
    Return _timeString($solNoonLocal, 3)
EndFunc   ;==>_calcSolNoon


Func _dayString($jd, $next, $flag)
    ; returns a string in the form DDMMMYYYY[ next] to display prev/next rise/set
    ; flag=2 for DD MMM, 3 for DD MM YYYY, 4 for DDMMYYYY next/prev
    If ($jd < 900000) Or ($jd > 2817000) Then
        Local $output = "error"
    Else
        Local $z = Floor($jd + 0.5)
        Local $f = ($jd + 0.5) - $z
        If $z < 2299161 Then
            $A = $z
        Else
            Local $alpha = Floor(($z - 1867216.25) / 36524.25)
            $A = $z + 1 + $alpha - Floor($alpha / 4)
        EndIf
        Local $B = $A + 1524
        Local $C = Floor(($B - 122.1) / 365.25)
        Local $D = Floor(365.25 * $C)
        Local $E = Floor(($B - $D) / 30.6001)
        Local $day = $B - $D - Floor(30.6001 * $E) + $f
        Local $month
        If $E < 14 Then
            $month = $E - 1
        Else
            $month = $E - 13
        EndIf
        Local $year
        If $month > 2 Then
            $year = $C - 4716
        Else
            $year = $C - 4715
        EndIf
        If $flag = 2 Then
            Local $output = _zeroPad($day, 2) & " " & $MonthAbrev[$month - 1]
        ElseIf $flag = 3 Then
            Local $output = _zeroPad($day, 2) & $MonthAbrev[$month - 1] & String($year)
        ElseIf $flag = 4 Then
            If $next Then
                Local $output = _zeroPad($day, 2) & $MonthAbrev[$month - 1] & String($year) & " next"
            Else
                Local $output = _zeroPad($day, 2) & $MonthAbrev[$month - 1] & String($year) & " prev"
            EndIf
        EndIf
    EndIf
    Return $output
EndFunc   ;==>_dayString


Func _TimeDateString($jd, $minutes)
    Local $output = _timeString($minutes, 2) & " " & _dayString($jd, 0, 2)
    Return $output
EndFunc   ;==>_TimeDateString


Func _timeString($minutes, $flag)
    ; timeString returns a zero-padded string (HH:MM:SS) given time in minutes
    ; flag=2 for HH:MM, 3 for HH:MM:SS
    If (($minutes >= 0) And ($minutes < 1440)) Then
        Local $floatHour = $minutes / 60.0
        Local $hour = Floor($floatHour)
        Local $floatMinute = 60.0 * ($floatHour - Floor($floatHour))
        Local $minute = Floor($floatMinute)
        Local $floatSec = 60.0 * ($floatMinute - Floor($floatMinute))
        Local $second = Floor($floatSec + 0.5)
        If ($second > 59) Then
            $second = 0
            $minute += 1
        EndIf
        If (($flag = 2) And ($second >= 30)) Then $minute += 1
        If ($minute > 59) Then
            $minute = 0
            $hour += 1
        EndIf
        Local $output = _zeroPad($hour, 2) & ":" & _zeroPad($minute, 2)
        If $flag > 2 Then $output = $output & ":" & _zeroPad($second, 2)
    Else
        $output = "error"
    EndIf
    Return $output
EndFunc   ;==>_timeString


Func _calcSunriseSetUTC($rise, $jd, $Latitude, $Longitude)
    Local $t = _CalcTimeJulianCent($jd)
    Local $eqTime = _calcEquationOfTime($t)
    Local $solarDec = _calcSunDeclination($t)
    Local $RiseSetAngle = _calcHourAngleSunrise($Latitude, $solarDec)
    If Not $rise Then $RiseSetAngle = -$RiseSetAngle
    Local $delta = $Longitude + _radToDeg($RiseSetAngle)
    Local $timeUTCRS = 720 - (4.0 * $delta) - $eqTime
    Return $timeUTCRS  ; in minutes
EndFunc   ;==>_calcSunriseSetUTC


Func _calcSunriseSet($rise, $jd, $Latitude, $Longitude, $TimeZone, $DST)
    ; rise = 1 for sunrise, 0 for sunset
    Local $timeUTC = _calcSunriseSetUTC($rise, $jd, $Latitude, $Longitude)
    Local $newTimeUTC = _calcSunriseSetUTC($rise, $jd + $timeUTC / 1440.0, $Latitude, $Longitude)
    If IsNumber($newTimeUTC) Then
        Local $timeLocal = $newTimeUTC + ($TimeZone * 60.0)
        If $rise Then
            Local $riseT = _CalcTimeJulianCent($jd + $newTimeUTC / 1440.0)
            Local $riseAz = _calcAzEl(0, $riseT, $timeLocal, $Latitude, $Longitude, $TimeZone)
        EndIf
        If $DST Then
            $timeLocal += 60.0
        Else
            $timeLocal += 0.0
        EndIf
        If (($timeLocal >= 0.0) And ($timeLocal < 1440.0)) Then
            Return _timeString($timeLocal, 2)
        Else
            Local $jday = $jd
            Local $increment
            If $timeLocal < 0 Then
                $increment = 1
            Else
                $increment = -1
            EndIf
            While (($timeLocal < 0.0) Or ($timeLocal >= 1440.0))
                $timeLocal += $increment * 1440.0
                $jday -= $increment
            WEnd
            Return _TimeDateString($jday, $timeLocal)
        EndIf
    Else
        Local $doy = _calcDoyFromJD($jd)
        If ((($Latitude > 66.4) And ($doy > 79) And ($doy < 267)) Or (($Latitude < -66.4) And (($doy < 83) Or ($doy > 263)))) Then
            If ($rise) Then
                $jdy = _calcJDofNextPrevRiseSet(0, $rise, $jd, $Latitude, $Longitude, $TimeZone, $DST)
            Else
                $jdy = _calcJDofNextPrevRiseSet(1, $rise, $jd, $Latitude, $Longitude, $TimeZone, $DST)
            EndIf
            ;Return _dayString ( $jdy, 0, 3 )
        Else
            If ($rise) Then
                $jdy = _calcJDofNextPrevRiseSet(1, $rise, $jd, $Latitude, $Longitude, $TimeZone, $DST)
            Else
                $jdy = _calcJDofNextPrevRiseSet(0, $rise, $jd, $Latitude, $Longitude, $TimeZone, $DST)
            EndIf
            ;Return _dayString ( $jdy, 0, 3 )
        EndIf
    EndIf
EndFunc   ;==>_calcSunriseSet


Func _calcJDofNextPrevRiseSet($next, $rise, $jd, $Latitude, $Longitude, $tz, $DST)
    Local $julianday = $jd
    Local $increment
    If $next Then
        $increment = 1.0
    Else
        $increment = -1.0
    EndIf
    Local $Time = _calcSunriseSetUTC($rise, $julianday, $Latitude, $Longitude)
    While Not IsNumber($Time)
        $julianday += $increment
        $Time = _calcSunriseSetUTC($rise, $julianday, $Latitude, $Longitude)
    WEnd
    Local $timeLocal
    If $DST Then
        $timeLocal = $Time + $tz * 60.0 + 60.0
    Else
        $timeLocal = $Time + $tz * 60.0
    EndIf
    While (($timeLocal < 0.0) Or ($timeLocal >= 1440.0))
        Local $incr
        If $timeLocal < 0 Then
            $incr = 1
        Else
            $incr = -1
        EndIf
        $timeLocal += ($incr * 1440.0)
        $julianday -= $incr
    WEnd
    Return $julianday
EndFunc   ;==>_calcJDofNextPrevRiseSet

;( $AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST )
Func _SolarCalculate($CalcLat, $CalcLong, $CalcMonth, $CalcDay, $CalcYear, $CalcTz, $CalcHour, $CalcMinute, $CalcSecond, $CalcDST, $CalcCurrent)
    Local $SolarCalculateResult[1] = [0]
    Local $jday = _getJD($CalcMonth, $CalcDay, $CalcYear)
    Local $tl = _getTimeLocal($CalcHour, $CalcMinute, $CalcSecond, $CalcDST)
    Local $tz
    If $CalcTz = "" Then
        $tz = Int($CalcLong / 15)
    Else
        $tz = $CalcTz
    EndIf
    Local $total = $jday + $tl / 1440.0 - $tz / 24.0
    Local $t = _CalcTimeJulianCent($total)
    Local $CalcSolNoon = _calcSolNoon($jday, $CalcLong, $tz, $CalcDST)    ; Returns String: "HH:MM:SS"
    ;MsgBox ( 0, "Solar Noon:", $CalcSolNoon )
    Local $rise = _calcSunriseSet(1, $jday, $CalcLat, $CalcLong, $tz, $CalcDST)    ; Returns String: "HH:MM" or "HH:MM DD Mon"
    ;MsgBox ( 0, "Sun Rise:", $rise )
    Local $set = _calcSunriseSet(0, $jday, $CalcLat, $CalcLong, $tz, $CalcDST)     ; Returns String: "HH:MM" or "HH:MM DD Mon"
    ;MsgBox ( 0, "Sun Set:", $set )
    _ArrayAdd($SolarCalculateResult, $rise)
    _ArrayAdd($SolarCalculateResult, $set)
    _ArrayAdd($SolarCalculateResult, $CalcSolNoon)
    $SolarCalculateResult[0] += 3
    If $CalcCurrent Then
        Local $CalcAzEl = _calcAzEl(1, $t, $tl, $CalcLat, $CalcLong, $tz)    ; Returns Array: 0 - Azimuth, 1 - Elevation, 2 - Illumination-Disposition
        _ArrayConcatenate($SolarCalculateResult, $CalcAzEl)
        $SolarCalculateResult[0] += 3
    EndIf
    Return $SolarCalculateResult
EndFunc   ;==>_SolarCalculate

Edited April 10, 2020 by SirAlonne
Clarification.

Melba23 · April 10, 2020

SirAlonne,

What does Google say for the dates in question? Does it produce similar times as well?

M23

SirAlonne · April 10, 2020

@Melba23
The results were the following:

UDF:
2020.03.30
Sunrise = 06:39 | Sunset = 19:31
------------------------------------------------
2020.03.31
Sunrise = 06:39 | Sunset = 19:31

Google:
2020.03.30
Sunrise = 06:39 | Sunset = 19:31
------------------------------------------------
2020.03.31
Sunrise = 06:39 | Sunset = 19:19

Sign In

Solar Position UDF - Debug

Recommended Posts

SirAlonne

Melba23

SirAlonne

Create an account or sign in to comment

Create an account

Sign in

Browse

AutoIt Resources

Release

Beta