こんにちは。
精度良好な地理的距離の短距離近似法として、Bowring's method for short lines (Wikipedia)1 を見つけました。
精度評価2を行うと、Gauss_mid_latitude(同じく短距離近似)よりも誤差は半球内では小さいです(下記のプロット内のbowring_short
)。
ソース
bowring_short.js
// Bowring, B. R. (1981). "The direct and inverse problems for short geodesic lines on the ellipsoid". Surveying and Mapping. 41 (2): 135–141.
function bowring_short(lat1, lon1, lat2, lon2){ // in degrees
"use strict";
const a = 6378137.0; // GRS80
const f = 1 / 298.257223563; // WGS84
const e2 = f * (2 - f); // first eccentricity squared
const ep2 = e2 / (1 - f) ** 2; // second eccentricity squared; epsilon
const degree = Math.PI / 180.0;
const sin = Math.sin;
const cos = Math.cos;
const sqrt = Math.sqrt;
const hypot = Math.hypot;
const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;
const londiffhalf = (lon2 - lon1) / 2 * degree;
lat1 = lat1 * degree;
lat2 = lat2 * degree;
const latdiffhalf = (lat2 - lat1) / 2;
const lat0 = lat1; // Bowring
const cos_lat0 = cos(lat0);
const cos_lat0_squared = cos_lat0 ** 2;
const C = sqrt(1 + ep2);
const B2 = 1 + ep2 * cos_lat0_squared;
const B = sqrt(B2);
const A = sqrt(1 + ep2 * cos_lat0_squared ** 2);
// radius, longitude, and latitude on the conformal sphere
const Rp = C / B2;
const lonpdiffhalf = A * londiffhalf;
const latpdiffhalf = latdiffhalf / B * (1 + (3 / 2 * ep2 / B2) * latdiffhalf * sin(2 * lat0 + 4 / 3 * latdiffhalf)); // D
const cos_latp0 = (B * cos_lat0 * cos(latpdiffhalf) - sin(lat0) * sin(latpdiffhalf)) / A;
// great circle formula
const E = cos(lonpdiffhalf) * sin(latpdiffhalf);
const F = sin(lonpdiffhalf) * cos_latp0;
return a * Rp * asin2(hypot(E, F)); // geodesic distance
}