地理的距離の簡素な短距離近似式の導出

こんにちは。
地理座標系（経緯度座標）で表わされた二点間 $\boldsymbol{p}_1, \ \boldsymbol{p}_2$の地理的距離の簡素な短距離近似式を導出し、精度評価を行いました。

### 近似式
便宜上、straight line 法と呼ぶことにしました。

```math
\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](https://en.wikipedia.org/wiki/Geographical_distance#Tunnel_distance) (Wikipedia) の計算を回転楕円体へ適用した形です。

```math
\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}
```

### 他の類似計算式
より精度向上した短距離近似式があります（下記）。

#### ガウスの中間緯度法

```math
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](https://en.wikipedia.org/wiki/Geographical_distance#Bowring's_method_for_short_lines) (Wikipedia)[^1] があります。

参考：[Bowring for short lines（精度良好な地理的距離の短距離近似法）](https://qiita.com/kkdd/items/3f56aed15301b24e8e5d)」、「[Bowring for short lines（精度良好な地理的距離の短距離近似法）をさらに精度改善（中間緯度を利用）](https://qiita.com/kkdd/items/f4ca002cc66637453f38)」

[^1]: [Bowring, B. R., "The direct and inverse problems for short geodesic lines on the ellipsoid". Surveying and Mapping](https://books.google.co.jp/books?id=eEI4AAAAIAAJ&pg=PA36&lpg=PA36&dq=Bowring+Surveying+and+Mapping+1981+degrees+radians&source=bl&ots=g_-TH3_hvQ&sig=ACfU3U01xtWe6_63kQoU_j2PKetjNk180Q&hl=ja&sa=X&ved=2ahUKEwjissjb9auHAxVoaPUHHfGoApwQ6AF6BAgVEAM#v=snippet&q=Bowring%20short%20geodesic%20ellipsoid&f=false) (1981), p.135.

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

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

```math
\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}
```

### 近似精度の比較評価[^2]
相対誤差最大値（縦軸）の距離依存性を評価しています（プロット内の```straightline```）。
![mini1736353739.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/54617/a1ff08aa-cb5a-195c-0f0d-11bd05a43726.jpeg)


[^2]: 参考： [二点間測地線距離（地球上）：各計算方法の誤差評価](https://qiita.com/kkdd/items/987e79d1bbf38fd84e71)

### source code
```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
  );
}
```