2地点からの方角で目標の位置を決める

緯度・経度によって与えられた2地点からの方角で目標点の緯度・経度を決める問題。ただし、方角というのは、その角度の測地線(大円)上に目標点がある、という意味である。「知恵袋」から拾った問題で興味をもったので、取り組んでみることにした。球面三角法でマニュアル化されるような問題であるかもしれないが、球面三角法は一から復習したくないので、3次元回転を用いることとした。

緯度・経度 (\lambda_\rm{A}, \mu_\rm{A}) の地点Aからの方角が東から反時計回りに \phi_\rm{A}、同様に (\lambda_\rm{B}, \mu_\rm{B}) の地点Bからの方角が \phi_\rm{B} である目標点の緯度・経度を求める。目標点は2つあるが、紛れのないようA,Bに近い方とする。図で \phi を目標の方角の反対側にとっているが、もともと2つの大円の交点座標が目的であったので、そうなっている。

まず、緯度・経度が (\lambda, \mu) の地点で方角 \phi の大円の方程式を求めることを考える。緯度・経度 (0, 0) における方角 0 の大円すなわち赤道を回転によって、目標とする大円に移す。

地球中心を原点とし、緯度・経度 (0, 0) に向かう方向を x-軸、極方向を z-軸とする。z-軸方向単位ベクトル \boldsymbol{e}_z を目標とする大円の法線方向に移動する。

z 軸まわりに角度 \theta 回転する回転行列は、

\sf{R}_z = \begin{pmatrix}\quad\cos\theta\qquad-\sin\theta\qquad 0\\\quad\sin\theta\qquad\quad\cos\theta \qquad 0\\\quad 0 \quad\quad\quad\quad 0 \quad\quad\quad 1\end{pmatrix}

と書ける。x, y 軸まわりも同様。ただし、回転の向きは右ねじの法則にしたがうものとする。

すると、目的の回転は

\sf{R}(\lambda, \mu, \phi) = \sf{R}_z(\mu) \sf{R}_y(-\lambda) \sf{R}_x(\phi)

となる。

この回転によって \boldsymbol{e}_z を目標の大円の法線単位ベクトル

\boldsymbol{e}_\zeta = \sf{R}(\lambda, \mu, \phi) \boldsymbol{e}_z

に移す。目標の大円は、この法線ベクトルに垂直だから、位置ベクトルを \boldsymbol{r} = (x, y, z) とするとき、

 \boldsymbol{r} \cdot \boldsymbol{e}_\zeta = 0

これが大円の方程式を表す。今回距離は目的でないので、地球半径は 1 として
x^2+y^2+z^2 = 1
を条件として加える。
2地点における \boldsymbol{e}_{\zeta \rm{A}}, \boldsymbol{e}_{\zeta \rm{B}} を求めれば、連立方程式

 \boldsymbol{r} \cdot \boldsymbol{e}_{\zeta \rm{A}} = 0
 \boldsymbol{r} \cdot \boldsymbol{e}_{\zeta \rm{B} }= 0

の解が目的とする地点の座標となる。

\lambda = \cos^{-1}\sqrt{x^2+y^2}
\mu = \tan^{-1}\displaystyle\frac{y}{x}

によって、その緯度、経度を得る。ただし、逆三角関数は値域に限定があるので、たとえば経度は上の値に 180° を足さねばならなかったりする。

以上の操作は、とても手計算でやろうとは思わないシロモノである。Mathcad2001 という年代物のソフトで行列計算をさせて正しく結果を出しているらしいことを確認した。購入時思い切った散財をした優秀な科学技術計算ソフトである。私はこのソフトを使うためにだけWindows XPマシンを1台残してある。現行の Mathcad は四半世紀を経て機能的にも進化をとげ、もはや価格において個人では到底手の出ないものとなっている。