Convert TWD97 to latitude, longitude

twd97.swift

// from https://gist.github.com/pingyen/1346895

import Foundation

extension Double {
    var ²: Double {
        get {
            return self * self
        }
    }
    
    var ³: Double {
        get {
            return self * self * self
        }
    }
    
    var ⁴: Double {
        get {
            return (self.²).²
        }
    }
    
    var ⁵: Double {
        get {
            return self * (self.²).²
        }
    }
    
    var ⁶: Double {
        get {
            return (self.³).²
        }
    }
}

func toLatitudeLongitude(x: Double, y: Double) -> (Double, Double) {
    let a: Double = 6378137.0
    let b: Double = 6356752.314245
    let lng0: Double = 121 * .pi / 180
    let k0: Double = 0.9999
    let dx: Double = 250000
    
    let dy: Double = 0.0
    let e: Double = sqrt(1.0 - (b * b) / (a * a))
    
    let x = x - dx
    let y = y - dy
    
    let M = y / k0

    let mu = M / (a * (1.0 - e.² / 4.0 - 3 * e.⁴ / 64.0 - e.⁶ / 256.0))
    
    let e1 = (1.0 - sqrt(1.0 - e.²)) / (1.0 + sqrt(1.0 - e.²))
    
    let J1 = 3 * e1 / 2 - 27 * e1.³ / 32
    let J2 = 21 * e1.² / 16 - 55 * e1.⁴ / 32
    let J3 = 151 * e1.³ / 96
    let J4 = 1097 * e1.⁴ / 512
    
    let fp = mu + J1 * sin(2 * mu) + J2 * sin(4 * mu) + J3 * sin(6 * mu) + J4 * sin(8 * mu)
    
    let e2 = (e * a / b).²
    let C1 = (e2 * cos(fp)).²
    let T1 = tan(fp).²
    let R1 = a * (1 - e.²) / pow(1 - e.² * sin(fp).², 3 / 2)
    let N1 = a / sqrt(1 - e.² * sin(fp).²)
    
    let D = x / (N1 * k0)
    
    let Q1 = N1 * tan(fp) / R1
    let Q2 = D.² / 2
    let Q3 = (5 + 3 * T1 + 10 * C1 - 4 * C1.² - 9 * e2) * D.⁴ / 24
    let Q4 = (61 + 90 * T1 + 298 * C1 + 45 * T1.² - 3 * C1.² - 252 * e2) * D.⁶ / 720
    
    let lat = fp - Q1 * (Q2 - Q3 + Q4)
    
    let Q5 = D
    let Q6 = (1 + 2 * T1 + C1) * D.³ / 6
    let Q7 = (5 - 2 * C1 + 28 * T1 - 3 * C1.² + 8 * e2 + 24 * T1.²) * D.⁵ / 120
    let lng = lng0 + (Q5 - Q6 + Q7) / cos(fp)
    
    return ((lat * 180) / .pi, (lng * 180) / .pi)
}