球面大円距離計算式の比較（浮動小数点計算精度）

こんにちは。
下記の球面大円距離計算式¹の浮動小数点計算精度を比較しました（近距離二点間、対蹠点間近傍の条件）。

• 球面三角法余弦公式による計算式（cosine_law）
• haversine 計算式
  　* 加えて、その堅牢化式（haversine_robust; atan2関数を利用）

短評

• cosine_law は、極端な近距離条件では、最大約 1 mm の誤差が生じますが、実用上は問題となることは少なそうです。
• haversine_robust は、haversine とほぼ差は現れず、出番はなさそうです。

計算例

$ node haversine.js 
赤道線上近距離二点間距離計算例 (m)
accurate: 1.1119492664455874
- error in haversine: 0
- error in haversine_robust: 2.220446049250313e-16
- error in cosine_law: -0.0007600928187647682
accurate: 0.0011119492664455873
- error in haversine: 0
- error in haversine_robust: 0
- error in cosine_law: -0.0011119492664455873

赤道線上対蹠点間近傍距離計算例 (m)
accurate: 20015085.684071306
- error in haversine: 0.004822958260774612
- error in haversine_robust: 0.004822958260774612
- error in cosine_law: 0.0007600933313369751
accurate: 20015086.79490862
- error in haversine: 0.0011119507253170013
- error in haversine_robust: 0.0011119507253170013
- error in cosine_law: 0.0011119507253170013

source code

haversine.js

'use strict'
const R_earth = 6371000;  // in m

function toRad(w) {  // in degrees
    return w * Math.PI / 180;
}

function haversine(lat1, lon1, lat2, lon2, R = R_earth) {  // in degrees
    const asin2 = w => w < 1 ? 2 * Math.asin(Math.sqrt(w)) : Math.PI;
    const φ1 = toRad(lat1), φ2 = toRad(lat2);
    const Δφhalf = (φ2 - φ1) / 2;
    const Δλhalf = toRad(lon2 - lon1) / 2;
    const y2 = Math.sin(Δφhalf) ** 2;
    const x2 = Math.sin(Δλhalf) ** 2 * Math.cos(φ1) * Math.cos(φ2);
    return R * asin2(x2 + y2);
}
//  not numerically robust extremely close to the antipodal point

function haversine_robust(lat1, lon1, lat2, lon2, R = R_earth) {  // in degrees
    const asin2 = w => w < 1 ? 2 * Math.atan2(Math.sqrt(w), Math.sqrt(1 - w)) : Math.PI;  // robust
    const φ1 = toRad(lat1), φ2 = toRad(lat2);
    const Δφhalf = (φ2 - φ1) / 2;
    const Δλhalf = toRad(lon2 - lon1) / 2;
    const y2 = Math.sin(Δφhalf) ** 2;
    const x2 = Math.sin(Δλhalf) ** 2 * Math.cos(φ1) * Math.cos(φ2);
    return R * asin2(x2 + y2);
}

function cosine_law_spherical_trigonometry(lat1, lon1, lat2, lon2, R = R_earth) {  // in degrees
    const Δλ = toRad(lon2 - lon1);
    const φ1 = toRad(lat1), φ2 = toRad(lat2);
    return R * Math.acos(Math.sin(φ1)*Math.sin(φ2)+Math.cos(φ1)*Math.cos(φ2)*Math.cos(Δλ));
}
//  not numerically robust for extremely close two points nor near the antipodal point

function distance_on_equator(lon1, lon2, R = R_earth) {
    return toRad(Math.abs(lon2 - lon1)) * R
}

const lat1 = 0.0, lon1 = 0.0, lat2 = 0.0;
let lon2, dist;
lon2 = 0.00001;
console.log("赤道線上近距離二点間距離計算例 (m)")
dist = distance_on_equator(lon1, lon2);
console.log("accurate:", dist)
console.log("error in haversine:", haversine(lat1, lon1, lat2, lon2) - dist);
console.log("error in haversine_robust:", haversine_robust(lat1, lon1, lat2, lon2) - dist);
console.log("error in cosine_law:", cosine_law_spherical_trigonometry(lat1, lon1, lat2, lon2) - dist);

lon2 = 180.0 - lon2;
console.log("赤道線上対蹠点間近傍距離計算例 (m)")
dist = distance_on_equator(lon1, lon2);
console.log("accurate:", dist)
console.log("error in haversine:", haversine(lat1, lon1, lat2, lon2) - dist);
console.log("error in haversine_robust:", haversine_robust(lat1, lon1, lat2, lon2) - dist);
console.log("error in cosine_law:", cosine_law_spherical_trigonometry(lat1, lon1, lat2, lon2) - dist);
console.log("")

lon2 = 0.00000001;
console.log("赤道線上近距離二点間距離計算例 (m)")
dist = distance_on_equator(lon1, lon2);
console.log("accurate:", dist)
console.log("error in haversine:", haversine(lat1, lon1, lat2, lon2) - dist);
console.log("error in haversine_robust:", haversine_robust(lat1, lon1, lat2, lon2) - dist);
console.log("error in cosine_law:", cosine_law_spherical_trigonometry(lat1, lon1, lat2, lon2) - dist);

lon2 = 180.0 - lon2;
console.log("赤道線上対蹠点間近傍距離計算例 (m)")
dist = distance_on_equator(lon1, lon2);
console.log("accurate:", dist)
console.log("error in haversine:", haversine(lat1, lon1, lat2, lon2) - dist);
console.log("error in haversine_robust:", haversine_robust(lat1, lon1, lat2, lon2) - dist);
console.log("error in cosine_law:", cosine_law_spherical_trigonometry(lat1, lon1, lat2, lon2) - dist);

1. 参考：https://ja.wikipedia.org/wiki/大円距離#数式 ↩