こんにちは。
地理座標系（経緯度座標）で表わされた二点間 $\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$です。

導出法

単純に、三次元デカルト座標系で直線距離（近似）を求めています。（新規性はないと思うのですが、同様の取り組みを見つけられませんでした）。Tunnel_distance (Wikipedia) の計算を回転楕円体へ適用した形です。

\begin{align}
& \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} \approx -M\left(\phi_\textrm{m}\right) \sin\phi_\textrm{m} \Delta \phi \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} \approx  2 N\left(\phi_\textrm{m}\right) \cos \phi _\textrm{m}  \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) \approx M\left(\phi_\textrm{m}\right) \cos\phi _\textrm{m} \Delta \phi ,\\
& {\left(D_\textrm{t} \left(\boldsymbol{p}_1 , \boldsymbol{p}_2 \right)\right)}^2 \triangleq (\Delta X)^2 + (\Delta Y)^2 + (\Delta Z)^2 \\
& \approx \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\right)^2 (\cos^2 \frac{\Delta \lambda}{2} + \cos^2 \phi_\textrm{m} \sin^2 \frac{\Delta \lambda}{2} ) \\
& \approx \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) \cos \frac{\Delta \lambda}{2} \Delta \phi\right)^2 .\\
\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

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

参考：Bowring for short lines（精度良好な地理的距離の短距離近似法）」、「Bowring for short lines（精度良好な地理的距離の短距離近似法）をさらに精度改善（中間緯度を利用）」

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

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

\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）。

ソースコード

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; // mid-latitude
  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
  );
}