Vincenty法は地球を回転楕円体として2点間の距離を近似する反復計算アルゴリズムである
このアルゴリズムには順解法と逆解法が存在する
始点座標と視点からの方位角、距離から終点座標を求めるアルゴリズム
2点の座標から方位角と距離を求めるアルゴリズム
逆解法に関しては本テーマとはずれているため説明はしない
本記事では以下のように定義している
$a=6378137.06$ 赤道半径
$f= \frac{1}{298.257223563}$ 扁平率
$b=(1-f)a$ 極半径
$\phi_{1},\phi_{2}$ 各点の緯度
$L_{1},L_{2}$ 各点の経度
$L=L_{2}-L_{1}$ 経度の差
$U_{1},U_{2}$ 補助球上の経度
$\lambda_{1},\lambda_{2}$ 補助球上の緯度
$\alpha_{1},\alpha_{2}$ 各点における方位角
$\alpha$ 赤道上での方位角
$s$ 2点間の距離
$ \sigma $ 補助球上の弧の長さ
まず初めに以下の計算を行う
$$
U_{1}=\arctan((1-f)\tan\phi_{1})\\
\sigma_{1} = \arctan \begin{pmatrix}\frac{\tan U_{1}}{\cos \alpha_{1}}\end{pmatrix}\\
\sin \alpha = \cos U_{1} \sin \alpha_{1}\\
\cos^2\alpha = 1-\sin^2\alpha\\
u^2 = \cos^2\alpha\frac{a^2-b^2}{b^2}\\
A=1+\frac{u^2}{16384} \lbrace 4096+u^2 \lbrack-768+u^2(320-175u^2)\rbrack\rbrace\\
B = \frac{u^2}{1024}\lbrace256+u^2\lbrack-128+u^2(74-47u^2)\rbrack\rbrace
$$
次に$\sigma=\frac{s}{bA}$で初期化し、以下の反復計算を収束するまで行う
$$
2\sigma_{m} = 2\sigma_{1}+\sigma\\
\varDelta \sigma=B\sin\sigma
\lbrace\cos(2\sigma_{m})+\frac{1}{4}B
\lbrack
\cos\sigma(-1+2\cos^2(2\sigma_{m}))-\frac{1}{6}B\cos(2\sigma_{m})(-3+4\sin^2\sigma)(-3+4\cos^2(2\sigma_{m}))
\rbrack
\rbrace\\
\sigma=\frac{s}{bA}+\varDelta\sigma
$$
$\sigma$が任意の精度まで収束したら以下の計算を行う
$$
\phi_{2}=\arctan
\begin{pmatrix}
\frac{
\sin U_{1}\cos\sigma+\cos U_{1}\sin\sigma\cos\alpha_{1}
}{
(1-f)\sqrt{\sin^2\alpha+(\sin U_{1}\sin\sigma-\cos U_{1}\cos\sigma\cos\alpha_{1})^2}
}
\end{pmatrix}\\
\lambda=\arctan
\begin{pmatrix}
\frac{
\sin\sigma\sin\alpha_{1}
}{
\cos U_{1}\cos\sigma-\sin U_{1}\sin\sigma\cos\alpha_{1}
}
\end{pmatrix}\\
C= \frac{f}{16}\cos^2\alpha
\lbrack
4+f(4-3\cos^2\alpha)
\rbrack\\
L = \lambda-(1-C)f\sin\alpha
\lbrace
\sigma+C\sin\sigma
\lbrack
\cos(2\sigma_{m})+C\cos\sigma(-1+2\cos^2(2\sigma_{m}))
\rbrack
\rbrace\\
L_{2} = L+L_{1}\\
\alpha_{2}=\arctan
\begin{pmatrix}
\frac{
\sin\alpha
}{
-\sin U_{1}\sin\sigma+\cos U_{1}\cos\sigma\cos\alpha_{1}
}
\end{pmatrix}
$$
自力で計算するのはなかなか大変ですが、「こうゆう方法があるよ」という程度で覚えておくと後々役に立つかもしれません
今後、このアルゴリズムを実装したプログラムも紹介するつもりです