こんにちは。
精度良好な地理的距離の短距離近似法として、Bowring's method for short lines (Wikipedia)¹ を見つけました。

精度評価²を行うと、Gauss_mid_latitude（同じく短距離近似）よりも誤差は表半球内では小さいです（下記のプロット内のbowring_short）。

ただし、裏半球では誤差の増大が顕著です。
続編でさらに精度改善しました：「Bowring for short lines（精度良好な地理的距離の短距離近似法）をさらに精度改善（中間緯度を利用）」

ソース

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
}

精度評価