SunのSwiftでの位置

Question

Swift3で太陽の位置を計算するためにthis solutionを実装しようとしています。私はこれを、10分ごとに23時50分まで踏み込んだ真夜中からの日を単に循環する別の関数でこれをラップします。

私は本当にRを理解していないし、私は完全に理解していないいくつかの詳細があります。特に、角かっこ付きのある種のif/clamp機能があるようです。私が混乱したときのPythonバージョンと比較して、私は最善を尽くしました。そうでなければ、唯一の違いはNSDateの使用によるものであり、これは先頭のコードの一部を単純化したものです。

戻ってくる値のいくつかは正しいように見え、結果をプロットするとカーブの基礎を見ることができます。しかし、1つのコール、つまり7AM、次に7:10のコールの結果は大きく異なります。

私は締め付けに何か間違っていると強く思っています。入力のマイナーな変更は、さまざまな方法でmod/truncedを取得し、出力をスイングさせます。しかし、私はそれを見つけることができません。これを理解している人は助けてもらえますか？

は、ここで私が取得しています出力のサンプルです：

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

[OK]を、OSXのためのRインタプリタをダウンロードした後、それがすべてで印刷を行うための複数の方法があることを発見し、何のデバッガを持っていないことを発見自分の警告、その他など、私が探していた問題が見つかりました。実際に値の1つを間違って締めていました。ここでは、C言語のような言語に簡単に変換でき、オリジナルより読みやすいSwift3のバージョンがあります。ターゲットプラットフォームの日付形式で機能する最初の2つの関数の独自のバージョンを提供する必要があります。 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 
    if h == 0 && m == 0 { 
     f = 0.0 
    } else if h == 0 { 
     f = m/60.0 
    } else if m == 0 { 
     f = h/24.0 
    } else { 
     f = (h/24.0) + (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) 
}