Sun 在 Swift 中的位置

Question

我正在尝试实现this solution 来计算 Swift3 中的太阳位置。然后，我将其包装在另一个函数中，该函数只是从午夜每 10 分钟一次到 23:50 循环一天。

我不太了解 R，并且我没有完全理解答案的一些细节，特别是似乎是某种带有方括号的 if/clamp 函数。我尽力了，当我感到困惑时与 Python 版本进行比较。否则唯一的区别是由于使用了NSDate，它简化了顶部的一些代码。

我得到的一些值看起来是正确的，当我绘制结果时我可以看到曲线的基础。但是，第一个电话（比如早上 7 点）和下一个电话（早上 7 点 10 分）的结果截然不同。

我强烈怀疑我在钳制方面做错了，输入中的微小变化会以不同的方式被修改/截断并摆动输出。但我看不出来。了解此算法的任何人都可以提供帮助吗？

这是我得到的输出示例：

2017-06-21 00:10:00 +0000 -16.0713262209521 31.7135341633943
2017-06-21 00:20:00 +0000 61.9971433936385 129.193513530349
2017-06-21 00:30:00 +0000 22.5263575559266 78.5445189561018
2017-06-21 00:40:00 +0000 29.5973897349096 275.081637736092
2017-06-21 00:50:00 +0000 41.9552795956374 262.989819486864

如您所见，它在迭代之间摇摆不定。地球不会那样转！我的代码如下，这个版本只是将结果发送到日志：

class func julianDayFromDate(_ date: Date) -> Double {
    let ti = date.timeIntervalSince1970
    return ((ti / 86400.0) + 2440587)
}

class func sunPath(lat: Double, lon: Double, size: CGSize) -> UIImage {
    var utzCal = Calendar(identifier: .gregorian)
    utzCal.timeZone = TimeZone(secondsFromGMT: 0)!
    let year = utzCal.component(.year, from: Date())
    let june = DateComponents(calendar: utzCal, year: year, month: 6, day: 21).date!

    // now we loop for every 10 minutes (2 degrees) and plot those points
    for time in stride(from:0, to:(24 * 60), by: 10) {
        let calcdate = june.addingTimeInterval(Double(time) * 60.0)
        let (alt, az) = sun(date: calcdate, lat: lat, lon: lon)
        print(calcdate, alt, az)
    }

class func sun(date: Date, lat: Double, lon: Double) -> (altitude: Double, azimuth: Double) {
    // these come in handy
    let twopi = Double.pi * 2
    let deg2rad = Double.pi / 180.0

    // latitude to radians
    let lat_radians = lat * deg2rad

    // the Astronomer's Almanac method used here is based on Epoch 2000, so we need to
    // convert the date into that format. We start by calculating "n", the number of
    // days since 1 January 2000
    let n = julianDayFromDate(date) - 2451545.0

    // it continues by calculating the position in ecliptic coordinates,
    // starting with the mean longitude of the sun in degrees, corrected for aberation
    var meanlong_degrees = 280.460 + (0.9856474 * n)
    meanlong_degrees = meanlong_degrees.truncatingRemainder(dividingBy: 360.0)

    // and the mean anomaly in degrees
    var meananomaly_degrees = 357.528 + (0.9856003 * n)
    meananomaly_degrees = meananomaly_degrees.truncatingRemainder(dividingBy: 360.0)
    let meananomaly_radians = meananomaly_degrees * deg2rad

    // and finally, the eliptic longitude in degrees
    var elipticlong_degrees = meanlong_degrees + (1.915 * sin(meananomaly_radians)) + (0.020 * sin(2 * meananomaly_radians))
    elipticlong_degrees = elipticlong_degrees.truncatingRemainder(dividingBy: 360.0)
    let elipticlong_radians = elipticlong_degrees * deg2rad

    // now we want to convert that to equatorial coordinates
    let obliquity_degrees = 23.439 - (0.0000004 * n)
    let obliquity_radians = obliquity_degrees * deg2rad

    // right ascention in radians
    let num = cos(obliquity_radians) * sin(elipticlong_radians)
    let den = cos(elipticlong_radians)
    var ra_radians = atan(num / den)
    ra_radians = ra_radians.truncatingRemainder(dividingBy: Double.pi)
    if den < 0 {
        ra_radians = ra_radians + Double.pi
    } else if num < 0 {
        ra_radians = ra_radians + twopi
    }
    // declination is simpler...
    let dec_radians = asin(sin(obliquity_radians) * sin(elipticlong_radians))

    // and from there, to local coordinates
    // start with the UTZ sidereal time
    let cal = Calendar.current
    let h = Double(cal.component(.hour, from: date))
    let m = Double(cal.component(.minute, from: date))
    let f: Double
    if h == 0 && m == 0 {
        f = 0.0
    } else if h == 0 {
        f = 60.0 / m
    } else if h == 0 {
        f = 24.0 / h
    } else {
        f = (24.0 / h) + (60.0 / m)
    }
    var utz_sidereal_time = 6.697375 + 0.0657098242 * n + f
    utz_sidereal_time = utz_sidereal_time.truncatingRemainder(dividingBy: 24.0)

    // then convert that to local sidereal time
    var localtime = utz_sidereal_time + lon / 15.0
    localtime = localtime.truncatingRemainder(dividingBy: 24.0)
    var localtime_radians = localtime * 15.0  * deg2rad
    localtime_radians = localtime.truncatingRemainder(dividingBy: Double.pi)

    // hour angle in radians
    var hourangle_radians =  localtime_radians - ra_radians
    hourangle_radians = hourangle_radians.truncatingRemainder(dividingBy: twopi)

    // get elevation in degrees
    let elevation_radians = (asin(sin(dec_radians) * sin(lat_radians) + cos(dec_radians) * cos(lat_radians) * cos(hourangle_radians)))
    let elevation_degrees = elevation_radians / deg2rad

    // and azimuth
    let azimuth_radians = asin( -cos(dec_radians) * sin(hourangle_radians) / cos(elevation_radians))

    // now clamp the output
    let azimuth_degrees: Double
    if (sin(dec_radians) - sin(elevation_radians) * sin(lat_radians) < 0) {
        azimuth_degrees = (Double.pi - azimuth_radians) / deg2rad
    } else if (sin(azimuth_radians) < 0) {
        azimuth_degrees = (azimuth_radians + twopi) / deg2rad
    } else {
        azimuth_degrees = azimuth_radians / deg2rad
    }

    return (elevation_degrees, azimuth_degrees)
}

Answer 1

好的，在为 OSX 下载了 R 解释器后，发现它没有调试器，发现有多种方法可以打印所有内容并带有自己的警告等，等等，我找到了我正在寻找的问题。它确实错误地钳制了其中一个值。这是一个有效的 Swift3 版本，应该很容易转换为任何类似 C 的语言，并且比原始版本更容易阅读。您必须提供您自己的前两个函数的版本，这些函数与您的目标平台的日期格式一起使用。 truncatingRemainer 是某人的疯狂想法，认为 Double 上不应该有 % 运算符，这是一个正常的 MOD。

// convinience method to return a unit-epoch data from a julian date
class func dateFromJulianDay(_ julianDay: Double) -> Date {
    let unixTime = (julianDay - 2440587) * 86400.0
    return Date(timeIntervalSince1970: unixTime)
}
class func julianDayFromDate(_ date: Date) -> Double {
    //==let JD = Integer(365.25 * (Y + 4716)) + Integer(30.6001 * (M +1)) +
    let ti = date.timeIntervalSince1970
    return ((ti / 86400.0) + 2440587.5)
}
// calculate the elevation and azimuth of the sun for a given date and location
class func sun(date: Date, lat: Double, lon: Double) -> (altitude: Double, azimuth: Double) {
    // these come in handy
    let twopi = Double.pi * 2
    let deg2rad = Double.pi / 180.0

    // latitude to radians
    let lat_radians = lat * deg2rad

    // the Astronomer's Almanac method used here is based on Epoch 2000, so we need to
    // convert the date into that format. We start by calculating "n", the number of
    // days since 1 January 2000. So if your date format is 1970-based, convert that
    // a pure julian date and pass that in. If your date is 2000-based, then
    // just let n = date
    let n = julianDayFromDate(date) - 2451545.0

    // it continues by calculating the position in ecliptic coordinates,
    // starting with the mean longitude of the sun in degrees, corrected for aberation
    var meanlong_degrees = 280.460 + (0.9856474 * n)
    meanlong_degrees = meanlong_degrees.truncatingRemainder(dividingBy: 360.0)

    // and the mean anomaly in degrees
    var meananomaly_degrees = 357.528 + (0.9856003 * n)
    meananomaly_degrees = meananomaly_degrees.truncatingRemainder(dividingBy: 360.0)
    let meananomaly_radians = meananomaly_degrees * deg2rad

    // and finally, the eliptic longitude in degrees
    var elipticlong_degrees = meanlong_degrees + (1.915 * sin(meananomaly_radians)) + (0.020 * sin(2 * meananomaly_radians))
    elipticlong_degrees = elipticlong_degrees.truncatingRemainder(dividingBy: 360.0)
    let elipticlong_radians = elipticlong_degrees * deg2rad

    // now we want to convert that to equatorial coordinates
    let obliquity_degrees = 23.439 - (0.0000004 * n)
    let obliquity_radians = obliquity_degrees * deg2rad

    // right ascention in radians
    let num = cos(obliquity_radians) * sin(elipticlong_radians)
    let den = cos(elipticlong_radians)
    var ra_radians = atan(num / den)
    ra_radians = ra_radians.truncatingRemainder(dividingBy: Double.pi)
    if den < 0 {
        ra_radians = ra_radians + Double.pi
    } else if num < 0 {
        ra_radians = ra_radians + twopi
    }
    // declination is simpler...
    let dec_radians = asin(sin(obliquity_radians) * sin(elipticlong_radians))

    // and from there, to local coordinates
    // start with the UTZ sidereal time, which is probably a lot easier in non-Swift languages
    var utzCal = Calendar(identifier: .gregorian)
    utzCal.timeZone = TimeZone(secondsFromGMT: 0)!
    let h = Double(utzCal.component(.hour, from: date))
    let m = Double(utzCal.component(.minute, from: date))
    let f: Double // universal time in hours and decimals (not days!)
    if h == 0 && m == 0 {
        f = 0.0
    } else if h == 0 {
        f = m / 60.0
    } else if m == 0 {
        f = h
    } else {
        f = h + (m / 60.0)
    }
    var utz_sidereal_time = 6.697375 + 0.0657098242 * n + f
    utz_sidereal_time = utz_sidereal_time.truncatingRemainder(dividingBy: 24.0)

    // then convert that to local sidereal time
    var localtime = utz_sidereal_time + lon / 15.0
    localtime = localtime.truncatingRemainder(dividingBy: 24.0)
    let localtime_radians = localtime * 15.0  * deg2rad

    // hour angle in radians
    var hourangle_radians =  localtime_radians - ra_radians
    hourangle_radians = hourangle_radians.truncatingRemainder(dividingBy: twopi)

    // get elevation in degrees
    let elevation_radians = (asin(sin(dec_radians) * sin(lat_radians) + cos(dec_radians) * cos(lat_radians) * cos(hourangle_radians)))
    let elevation_degrees = elevation_radians / deg2rad

    // and azimuth
    let azimuth_radians = asin( -cos(dec_radians) * sin(hourangle_radians) / cos(elevation_radians))

    // now clamp the output
    let azimuth_degrees: Double
    if (sin(dec_radians) - sin(elevation_radians) * sin(lat_radians) < 0) {
        azimuth_degrees = (Double.pi - azimuth_radians) / deg2rad
    } else if (sin(azimuth_radians) < 0) {
        azimuth_degrees = (azimuth_radians + twopi) / deg2rad
    } else {
        azimuth_degrees = azimuth_radians / deg2rad
    }

    // all done!
    return (elevation_degrees, azimuth_degrees)
}