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

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

ソース

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 latdiff = (lat2 - lat1);
  const cos_lat1 = cos(lat1);
  const cos_lat1_squared = cos_lat1 ** 2;

  const C = sqrt(1 + ep2);
  const B = sqrt(1 + ep2 * cos_lat1_squared);
  const A = sqrt(1 + ep2 * cos_lat1_squared ** 2);

 // radius, longitude, and latitude on the conformally mapped sphere
  const Rp = C / B ** 2;
  const lonpdiffhalf = A * londiffhalf;
  const latpdiffhalf = latdiff / B / 2 * (1 + (3 * ep2 / 4 / B ** 2) * latdiff * sin(2 * lat1 + 2 / 3 * latdiff));
  const cos_latp = (B * cos_lat1 * cos(latpdiffhalf) - sin(lat1) * sin(latpdiffhalf)) / A;

  const E = cos(lonpdiffhalf) * sin(latpdiffhalf);
  const F = sin(lonpdiffhalf) * cos_latp;
  return a * Rp * asin2(hypot(E, F));  // geodesic distance
}

精度評価