こんにちは。
地理座標系（経緯度座標）で表わされた二点間 $\boldsymbol{p}_1, \ \boldsymbol{p}_2$の地理的距離の簡素な短距離近似式を導出しました（便宜上、straight line 法と呼ぶことにしました）。（下記で精度評価もしています。）

近似式

\begin{align}
& D \left(\boldsymbol{p}_1 , \boldsymbol{p}_2 \right) \approx \sqrt{\left(2 N\left(\phi_\textrm{m}\right) \cos\phi_\textrm{m} \sin \frac{\Delta \lambda}{2}\right)^2 + \left(M\left(\phi_\textrm{m}\right) \, \Delta \phi \, \cos \frac{\Delta \lambda}{2} \right)^2 } ,\\
& \boldsymbol{p}_1 \triangleq (\frac{\Delta \lambda}{2},\ \phi_\textrm{m} + \frac{\Delta \phi}{2},\ 0), \ \boldsymbol{p}_2 \triangleq (- \frac{\Delta \lambda}{2},\ \phi_\textrm{m} - \frac{\Delta \phi}{2},\ 0), \\
& N(\phi) = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi}}, \ M(\phi) = N(\phi) \frac{1 - e^2}{1 - e^2 \sin^2 \phi} \\
& e^2 = f \left(2 - f\right), \ f = 1 / 298.257223563 , \ a = 6378137.0 \, \text{m}
\end{align}

導出に用いた短距離近似条件は、$|\cos\phi_\textrm{m} \sin \frac{\Delta \lambda}{2}| \ll 1$ かつ $|\frac{\Delta \phi}{2}| \ll 1$ （および $\Delta \phi \cos\phi_\textrm{m} \sin \frac{\Delta \lambda}{2} \approx 0$）です。

これの誤差評価は、$| D_\text{approx} - D_\text{true}| \lessapprox \frac{D}{24} \left(\frac{D}{a}\right)^2$です。

導出法

単純に、三次元デカルト座標系で直線距離を求めています。（新規性はないと思うのですが、同様の取り組みを見つけられませんでした）。

\begin{align}
& D \left(\boldsymbol{p}_1 , \boldsymbol{p}_2 \right) \approx \sqrt{(\Delta X)^2 + (\Delta Y)^2 + (\Delta Z)^2} \\
& \Delta X = \left( N(\phi_\textrm{m} + \frac{\Delta \phi}{2})\cos(\phi_\textrm{m} + \frac{\Delta \phi}{2}) - N(\phi_\textrm{m} - \frac{\Delta \phi}{2})\cos(\phi_\textrm{m} - \frac{\Delta \phi}{2}) \right) \cos\frac{\Delta \lambda}{2},\\
& \Delta Y = \left( N(\phi_\textrm{m} + \frac{\Delta \phi}{2})\cos(\phi_\textrm{m} + \frac{\Delta \phi}{2}) + N(\phi_\textrm{m} - \frac{\Delta \phi}{2})\cos(\phi_\textrm{m} - \frac{\Delta \phi}{2}) \right) \sin\frac{\Delta \lambda}{2},\\
& \Delta Z = (1-e^2) \left( N(\phi _\textrm{m} + \frac{\Delta \phi}{2})\sin(\phi _\textrm{m} + \frac{\Delta \phi}{2}) - N(\phi _\textrm{m} - \frac{\Delta \phi}{2})\sin(\phi _\textrm{m} - \frac{\Delta \phi}{2}) \right) ,\\
\end{align}

他の類似計算式

なおより精度向上した短距離近似式として、「ガウスの平均緯度法」があります：

D \left(\boldsymbol{p}_1 , \boldsymbol{p}_2 \right) \approx 2 N\left(\phi_\textrm{m}\right) \arcsin \sqrt{\left(\sin \frac{\Delta \lambda}{2} \cos\phi_\textrm{m} \right)^2 + \left(\cos \frac{\Delta \lambda}{2} \sin \frac{M\left(\phi_\textrm{m}\right) \Delta \phi}{2 N\left(\phi_\textrm{m}\right)} \right)^2}

これの誤差評価は、$| D_\text{approx} - D_\text{true}| \lessapprox \frac{D}{400} \left(\frac{D}{a}\right)^2$です。

さらにより精度良好な短距離近似法として、Bowring's method for short lines (Wikipedia)¹ があります。

簡素さを優先した短距離近似式

ところで精度を追求するよりも、簡素さを優先した形を欲しい場合には、高緯度（特に極付近）を適用除外することが許されるならば、$|\Delta \lambda| \ll 1$としてよく、より簡素な短距離近似計算式へ帰着されます（広く使われ、平面法などと呼ばれています）：

\begin{align}
& D \left(\boldsymbol{p}_1 , \boldsymbol{p}_2 \right) \approx \sqrt{\left(N\left(\phi_\textrm{m}\right) \cos\phi_\textrm{m} \Delta \lambda\right)^2 + \left(M\left(\phi_\textrm{m}\right) \Delta \phi\right)^2 } 
\end{align}

近似精度の比較評価²

相対誤差最大値（縦軸）の距離依存性を評価しています。

ソース

straightline.js

function straightline(lat1, lon1, lat2, lon2){  // in degrees
  "use strict";
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const lat = (lat1 + lat2) / 2 * degree; // middle
  const sinlat = sin(lat);
  const coslat = cos(lat);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio	

  return 2 * n * a * hypot(
    sin(londiffhalf) * coslat,
    cos(londiffhalf) * latdiffhalf * m_by_n
  );
}