こんにちは。
精度良好な地理的距離の短距離近似法である Bowring's method for short lines (Wikipedia)¹ を、さらに精度改善しました。

用いた計算方法は、"Bowring for short lines"(Karney; geographiclib) をそっくり踏襲し、中間緯度（mid-latitude）を用いました。

精度評価²を行うと、

「パラメトリック緯度＋中間緯度」を用いた計算は、短距離近似法としては確かに最も精度良好のようで（下記のプロット内のbowring_short_beta_；100kmでは10ppb）、
「地理緯度＋中間緯度」を用いる計算はこれに次ぐようです（bowring_short_）。

精度評価

source code

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
}

// by Karney
// https://geographiclib.sourceforge.io/C++/doc/geodesic.html#bowring
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 tan = Math.tan;
  const atan = Math.atan;
  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 latm = (lat1 + lat2) / 2;  // mid-latitude by Karney
  const cosm2 = cos(latm) ** 2;

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

// radius, longitude, and latitude on the conformal sphere
  const Rp = C / B2;
  const lonpdiffhalf = A * londiffhalf;
  
  const dlat = lat1 - latm;  // = -(lat2 - lat)
  const phipm = atan(tan(latm)/B) + dlat*dlat * 3*ep2/(4*B*B2) * sin(2*latm) * cos(2/3*dlat);
  const phipdiffhalf = dlat/B * (1 + (3*ep2*dlat)/(4*B2) * cos(2*latm) * sin(2/3*dlat));

// great circle formula
  const E = cos(lonpdiffhalf) * sin(phipdiffhalf);
  const F = sin(lonpdiffhalf) * cos(phipm);
  return a * Rp * asin2(hypot(E, F));  // geodesic distance
}

// by Karney
// https://geographiclib.sourceforge.io/C++/doc/geodesic.html#bowring
function bowring_short_beta_(lat1, lon1, lat2, lon2){  // in degrees
  "use strict";
  const a = 6378137.0;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const f1 = 1 - f;
  const e2 = f * (2 - f);  // first eccentricity squared
  const ep2 = e2 / (f1*f1);  // second eccentricity squared; epsilon

  const degree = Math.PI / 180.0;
  const sin = Math.sin;
  const cos = Math.cos;
  const tan = Math.tan;
  const atan = Math.atan;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;
  const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;

  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const bet1 = atan(tan(lat1 * degree) * f1);  // reduced latitude
  const bet2 = atan(tan(lat2 * degree) * f1);
  const betm = (bet1 + bet2)/2;  // mid-latitude by Karney
  const sbetm = sin(betm);
  const cbetm = cos(betm);
  const phim = atan(tan(betm) / f1);
  const cosm2 = cos(phim) ** 2;
  const C = sqrt(1 + ep2);
  const B2 = 1 + ep2 * cosm2;
  const B = sqrt(B2);
  const A = sqrt(1 + ep2 * cosm2 ** 2);

// radius, longitude, and latitude on the conformal sphere
  const Rp = C / B2;
  const lonpdiffhalf = A * londiffhalf;
  const dnm2 = 1 + ep2*sbetm*sbetm;
  const dnm = sqrt(dnm2);
  const dbet = (bet1-betm)/dnm;  // = -(bet2-betm)/dnm
  const phipm = atan(tan(phim)/B) + dbet * dbet * ep2/2 * cbetm*sbetm/dnm;
  const phipdiffhalf = dbet * (1 + dbet * dbet * ep2/6 * (cbetm*cbetm/dnm2 - sbetm*sbetm) );

// great circle formula
  const E = cos(lonpdiffhalf) * sin(phipdiffhalf);
  const F = sin(lonpdiffhalf) * cos(phipm);
  return a * Rp * asin2(hypot(E, F));  // geodesic distance
}