 
         /* Compute the local sidereal time of this moment */
         let siderealTime = (GMST0(day: d) + Degree(180.0) + positiveEastwardLongitude).reduced
-        
+                
         /* Compute time when Sun is at south - in hours UT */
         let sun = Sun(julianDay: StandardEpoch_J2000_0 + JulianDay(d))
         let sunRa = sun.equatorialCoordinates.rightAscension.inDegrees.reduced0

            return (nil, nil, .alwaysBelowAltitude)
        }

        let jd_start = self.julianDay.midnight + 0.5
        var jd_rise = jd_start
        var jd_set  = jd_start
        let accuracy = 1.0.seconds.inDays
        var jd_diff  = 10*accuracy
        
        let radiusVector = self.radiusVector
        let UT_sun_crossing_meridian = find_time_Sun_at_South(forJulianDay: julianDay, coordinates: coordinates)
        
        var security = 0
        while abs(jd_diff) > accuracy && security < 10 {
            security += 1

            let arc1 = diurnalArc(forJulianDay: jd_rise, sunAltitude: sunAltitude, coordinates: coordinates, radiusVector: radiusVector)
            let arc2 = diurnalArc(forJulianDay: jd_set, sunAltitude: sunAltitude, coordinates: coordinates, radiusVector: radiusVector)
            
            if arc1.error != nil || arc2.error != nil {
                security = 10
                break
            }
            
            let new_jd_rise = Calendar.gregorianGMT.date(bySettingHour: UT_sun_crossing_meridian - arc1.value!, of: jd_rise.date).julianDay
            jd_set = Calendar.gregorianGMT.date(bySettingHour: UT_sun_crossing_meridian + arc2.value!, of: jd_set.date).julianDay

            jd_diff = new_jd_rise-jd_rise
            jd_rise = new_jd_rise

    }
}

func find_time_Sun_at_South(forJulianDay jd: JulianDay, coordinates: GeographicCoordinates) -> Hour {
    // Somehow, there is a 0.5 day shift (midnight instead of noon) compared to explanations given by Paul Schylter.
    // See http://www.stjarnhimlen.se/comp/riset.html
    // But the values are correct like this.
    let localMidnight = jd.localMidnight(longitude: coordinates.longitude)
    let siderealTime = localMidnight.meanLocalSiderealTime(longitude: coordinates.longitude)
    let sun = Sun(julianDay: localMidnight)
    let sunRa = sun.equatorialCoordinates.rightAscension.inDegrees.reduced0
    return Hour((siderealTime.inDegrees-sunRa).reduced0.value/15.04107)    
}
