Swift4 经纬度计算日出日落时间

import UIKit

class PGSun: NSObject {
    struct Sunriseset {
        var sunrise = "07:00" // 日出
        var sunset = "18:30" // 日落
        init(sunrise:String, sunset:String) {
            self.sunrise = sunrise
            self.sunset = sunset
        }
    }
    var start: Double = 0.0
    var sRA: Double = 0.0
    var sdec: Double = 0.0
    var sr: Double = 0.0
    var lon: Double = 0.0
    var end: Double = 0.0
    
    private let Inv360: Double = 1.0 / 360.0
    private let Degrad: Double = .pi / 180.0
    private let Radge: Double = 180.0 / .pi
    
    /// 根据经纬度计算日出和日落时间  -误差2分钟
    ///
    /// - Parameters:
    /// - Returns: Sunriseset
    class func sun(longitude:Double, latitude:Double)->Sunriseset{
        let sun = PGSun()
        return sun.getSunTimeAtDate(d: Date(), longitude: longitude, latitude: latitude)
    }
    
    func getSunTimeAtDate(d:Date , longitude:Double ,latitude:Double)->Sunriseset{
        let xcts: Int = days(since_2000_Jan_0: d)
        let dic = getSunTime(xcts, lon: longitude, lat: latitude)
        return dic
    }
    
    func getSunTime(_ day: Int, lon longitude: Double, lat latitude: Double) -> Sunriseset {
        _ = self.sunRiset(day, long:longitude, lat:latitude ,altit: -35.0 / 60.0, limb:1, trise:self.start, tset:self.end)
        let sunrise = toLocalTime(self.start)
        let sunset = toLocalTime(self.end)
        let sun = Sunriseset(sunrise: sunrise, sunset: sunset)
        return sun
    }
    func toLocalTime(_ utTime:Double)-> String {
        var hour = Int(floor(utTime))
        var temp = utTime - Double(hour)
        hour += 8;
        temp = temp * 60;
        let minute = Int(floor(temp))
        return String(format: "%02d:%02d", hour,minute)
    }
    
    func days(since_2000_Jan_0 d: Date) -> Int {
        //距离2000-01-01的天数
        let dateStr = "2000-01-01"
        let dateFormatter = DateFormatter()
        dateFormatter.dateFormat = "yyyy-MM-dd"
        let date: Date? = dateFormatter.date(from: dateStr)
        
        var time: TimeInterval? = nil
        if let aDate = date {
            time = d.timeIntervalSince(aDate)
        }
        let days = (Int(time ?? 0)) / (3600 * 24)
        return days
    }
    func sunRiset(_ day:Int, long longitude:Double, lat:Double ,altit:Double, limb upper_limb:Int, trise:Double, tset:Double)->Int{
        var altit  = altit
        var d = 0.0/* Days since 2000 Jan 0.0 (negative before) */
        //以历元2000.0起算的日数。
        var sradius = 0.0   /* Sun‘s apparent radius */
        //太阳视半径，约16分（受日地距离、大气折射等诸多影响）
        
        var  t = 0.0         /* Diurnal arc */
        //周日弧，太阳一天在天上走过的弧长。
        
        var  tsouth = 0.0    /* Time when Sun is at south */
        var sidtime = 0.0    /* Local sidereal time */
        //当地恒星时，即地球的真实自转周期。比平均太阳日（日常时间）长3分56秒。
        
        var rc = 0; /* Return cde from function - usually 0 */
        /* Compute d of 12h local mean solar time */
        d = Double(day)/* Days_since_2000_Jan_0(date)*/ + 0.5 - longitude / 360.0;
        //计算观测地当日中午时刻对应2000.0起算的日数。
        /* Compute local sideral time of this moment */
        sidtime = self.revolution(self.GMST0(d) + 180.0 + longitude)
        //计算同时刻的当地恒星时（以角度为单位）。以格林尼治为基准，用经度差校正。
        
        /* Compute Sun‘s RA + Decl at this moment */
        self.Sun_RA_dec(d, ra:self.sRA, dec:self.sdec, r:self.sr);
        //    Sun_RA_dec(d, sRA,sdec,sr);
        //计算同时刻太阳赤经赤纬。
        
        /* Compute time when Sun is at south - in hours UT */
        tsouth = 12.0 - self.rev180(sidtime - self.sRA) / 15.0;
        //    tsouth = 12.0 - Rev180(sidtime - sRA) / 15.0;
        //计算太阳日的正午时刻，以世界时（格林尼治平太阳时）的小时计。
        
        /* Compute the Sun‘s apparent radius, degrees */
        sradius = 0.2666 / self.sr;
        //太阳视半径。0.2666是一天文单位处的太阳视半径（角度）。
        
        /* Do correction to upper limb, if necessary */
        if upper_limb != 0 {
            altit -= sradius;
            //如果要用上边缘，就要扣除一个视半径。
        }
        /* Compute the diurnal arc that the Sun traverses to reach */
        //计算周日弧。直接利用球面三角公式。如果碰到极昼极夜问题，同前处理。
        /* the specified altitide altit: */
        
        var cost:Double = 0
        cost = (self.sind(altit) - self.sind(lat)*self.sind(self.sdec))/(self.cosd(lat) * self.cosd(self.sdec));
        //    cost = (Sind(altit) - Sind(lat) * Sind(sdec)) /
        //    (Cosd(lat) * Cosd(sdec));
        
        if (cost >= 1.0)
        {
            rc = -1;
            t = 0.0;
        }
        else
        {
            if (cost <= -1.0)
            {
                rc = +1;
                t = 12.0;      /* Sun always above altit */
            }
            else{
                //            t = Acosd(cost) / 15.0;   /* The diurnal arc, hours */
                t = self.acosd(cost) / 15.0;
            }
        }
        
        /* Store rise and set times - in hours UT */
        self.start = tsouth - t;
        self.end = tsouth + t;
        return rc;
    }
    
    func Sun_RA_dec(_ d: Double, ra RA: Double, dec: Double, r: Double) {
        var obl_ecl: Double = 0
        var x: Double = 0
        var y: Double = 0
        var z: Double = 0
        sunpos(d, lon: lon, r: r)
        //计算太阳的黄道坐标。
        x = sr * cosd(lon)
        y = sr * sind(lon)
        //计算太阳的直角坐标。
        obl_ecl = 23.4393 - 3.563E-7 * d
        //黄赤交角，同前。
        z = y * sind(obl_ecl)
        y = y * cosd(obl_ecl)
        //把太阳的黄道坐标转换成赤道坐标（暂改用直角坐标）。
        sRA = atan2d(y, x: x)
        sdec = atan2d(z, x: sqrt(x * x + y * y))
        //最后转成赤道坐标。显然太阳的位置是由黄道坐标方便地直接确定的，但必须转换到赤
        //道坐标里才能结合地球的自转确定我们需要的白昼长度。
        
    }
    func sunpos(_ d: Double, lon: Double, r: Double) {
        var lon = lon
        var M: Double = 0
        var w: Double = 0
        var e: Double = 0
        var E: Double = 0
        var x: Double = 0
        var y: Double = 0
        var v: Double = 0 //真近点角，太阳在任意时刻的真实近点角。
        
        M = revolution(356.0470 + 0.9856002585 * d)
        //自变量的组成：2000.0时刻太阳黄经为356.0470度,此后每天约推进一度（360度/365天
        w = 282.9404 + 4.70935E-5 * d //近日点的平均黄经。
        
        e = 0.016709 - 1.151E-9 * d //地球公转椭圆轨道离心率的时间演化。以上公式和黄赤交角公式一样，不必深究。
        
        E = M + e * Double(Radge) * sind(M) * (1.0 + e * cosd(M))
        x = cosd(E) - e
        y = sqrt(1.0 - e * e) * sind(E)
        sr = sqrt(x * x + y * y)
        v = atan2d(y, x: x)
        
        lon = v + w
        self.lon = lon
        if lon >= 360.0 {
            lon -= 360.0
            self.lon = lon
        }
    }
    func revolution(_ x: Double) -> Double {
        return x - 360.0 * floor(x * Double(Inv360))
    }
    func rev180(_ x: Double) -> Double {
        return x - 360.0 * floor(x * Double(Inv360) + 0.5)
    }
    
    func GMST0(_ d: Double) -> Double {
        var sidtim0: Double
        sidtim0 = revolution((180.0 + 356.0470 + 282.9404) + (0.9856002585 + 4.70935E-5) * d)
        return sidtim0
    }
    func sind(_ x: Double) -> Double {
        return sin(x * Double(Degrad))
    }
    func cosd(_ x: Double) -> Double {
        return cos(x * Double(Degrad))
    }
    func acosd(_ x: Double) -> Double {
        return Double(Radge * acos(x))
    }
    func atan2d(_ y: Double, x: Double) -> Double {
        return Double(Radge * atan2(y, x))
    }

}