ヒュベニの距離計算式による2地点の距離(m)を計算するコード。測地系にはWGS84を使用

hubeny-distance.js

/**
 * ヒュベニの距離計算式を使用して、2点({latitude1,longitude1}, {latitude2,longitude2})間の
 * 距離を求める。測地系にはWGS84を使用。
 *
 * 参考：
 * - ヒュベニの距離計算式の簡略式( http://www.kashmir3d.com/kash/manual/std_siki.htm )
 * - 正式な公式は http://www.amano-tec.com/apps/paceruler.html を参照
 */
function calculateDistance(latitude1, longitude1, latitude2, longitude2) {
    // 先に計算しておいた定数
    let e2 = 0.00669437999019758;   // WGS84における「離心率e」の2乗
    let Rx = 6378137.0;             // WGS84における「赤道半径Rx」
    let m_numer = 6335439.32729246; // WGS84における「子午線曲率半径M」の分子(Rx(1-e^2))

    function deg2rad(deg) {
        return deg * Math.PI / 180.0;
    }

    let rad_lat1 = deg2rad(latitude1);
    let rad_lon1 = deg2rad(longitude1);
    let rad_lat2 = deg2rad(latitude2);
    let rad_lon2 = deg2rad(longitude2);

    let dp = rad_lon1 - rad_lon2;       // 2点の緯度差
    let dr = rad_lat1 - rad_lat2;       // 2点の経度差
    let p = (rad_lon1 + rad_lon2) / 2.0;// 2点の平均緯度

    let w = Math.sqrt(1.0 - e2 * Math.pow(Math.sin(p), 2));
    let m = m_numer / Math.pow(w, 3);   // 子午線曲率半径
    let n = Rx / w;                     // 卯酉(ぼうゆう)線曲率半径

    // 2点間の距離(単位m)
    let d = Math.sqrt(Math.pow((m * d_lon), 2)
                      + Math.pow((n * Math.cos(p) * dr), 2));
    return d;
}