楕円積分を用いた楕円弧長の計算　その3

1. はじめに
   　前回は、弧の存在範囲を、第I象限から楕円の右半分に拡げた。その際楕円積分の$\varphi$の範囲を$[0,\pi/2]$に制限していた為、$\theta$の範囲に応じてやや複雑な場合分けと計算式が必要になった。今回は、弧の存在範囲を楕円全体に拡げ、同時に楕円積分と変換関数の拡張により、場合分けの煩わしさの軽減も図る。

2. 第2種楕円積分の拡張

積分の基本性質により、$\varphi$が負の時、$E_0(\varphi)=-E_0(-\varphi)$が、自動的に成立する。

拡張形と表現しているが、これは第2種不完全楕円積分の本来持っている基本的性質を、今回必要になる範囲について、$E_0(\varphi,k)$で表したものに過ぎない。$E(\varphi,k)$のもっと広い性質については、 DLMFのサイト などを参照されたい。(特に、19.2.10に注目。)

3. 変換関数の拡張
   $\theta$の範囲を$-\pi \leq \theta \leq \pi$に拡げ、これに合わせて変換関数を拡張する。
   $\hspace{10pt}$ $th2u(\theta) = \begin{eqnarray} \left\{ \begin{array}{l} \tan^{-1}((b/a)/\tan \theta) &( 0 \lt \theta \lt \pi,\, \theta \neq \pi/2 ) \\ \tan^{-1}((b/a)/\tan \theta)+\pi &( -\pi \lt \theta \lt 0,\, \theta \neq -\pi/2 ) \\ 3\pi/2 &( \theta = -\pi ) \\ \pi &( \theta = -\pi/2 ) \\ \pi/2 &( \theta = 0 ) \\ 0 &( \theta = \pi/2 ) \\ -\pi/2 &( \theta = \pi ) \end{array} \right. \end{eqnarray} $
   $th2u(\theta)$は、$[-\pi,\pi]$の範囲の$\theta$に対して$[-\pi/2,3\pi/2]$の範囲を(単調減少に)返す。
   $\hspace{10pt}$ $th2v(\theta) = \begin{eqnarray} \left\{ \begin{array}{l} \tan^{-1}((a/b)\tan \theta) &( -\pi/2 \lt \theta \lt \pi/2 ) \\ \tan^{-1}((a/b)\tan \theta)-\pi &( -\pi \lt \theta \lt -\pi/2) \\ \tan^{-1}((a/b)\tan \theta)+\pi &( \pi/2 \lt \theta \lt \pi ) \\ \theta &( \theta = -\pi, -\pi/2, \pi/2, \pi ) \end{array} \right. \end{eqnarray} $
   $th2v(\theta)$は、$[-\pi,\pi]$の範囲の$\theta$に対して$[-\pi,\pi]$の範囲を(単調増加に)返す。
   変換関数がだいぶ複雑になったので、結果を確認するための逆変換関数を導入しておく。
   $\hspace{10pt}$ $u2th(u) = atan2(b\cos u,\,a\sin u)$
   $\hspace{10pt}$ $v2th(v) = atan2(b\sin v,\,a\cos v)$
   ここで用いたatan2は、殆どすべてのプログラミング言語でこの名前で提供されている、2変数のarc-tangentで、普通のarc-tangentで表すと次のようになる。
   $\hspace{10pt}$ $ atan2(y,x) = \begin{eqnarray} \left\{ \begin{array}{l} tan^{-1}(y/x) &(x \gt 0) \\ tan^{-1}(y/x) + \pi &(x\lt 0,\,y \geq 0 ) \\ tan^{-1}(y/x) - \pi &(x\lt 0,\,y \lt 0) \\ \pi/2 &(x = 0,\,y \gt 0) \\ -\pi/2 &(x = 0,\,y \lt 0)　\\ 0 &(x = 0,\,y = 0) \end{array} \right. \end{eqnarray} $

4. 楕円の弧長
   拡張された楕円積分と変換関数を使う事で、$\theta$の範囲による場合分けが不要になっている。
   $\theta$を動径(楕円上の点$(x,y)$と原点を結ぶ線分)が$x$軸となす角度とする。（反時計方向を正とする。）
   この時、$\theta$の$[\theta_1,\theta_2](-\pi \leq \theta_1 \leq \theta_2 \leq \pi)$の範囲の弧長は、次のように求められる。
   A) $0 \lt b \leq a$の場合
   $\hspace{10pt}$ $e = \sqrt{1-(b/a)^2}$
   $\hspace{10pt}$ $ u_1 = th2u(\theta_2)$
   $\hspace{10pt}$ $ u_2 = th2u(\theta_1)$
   とすると、弧長$L$は
   $\hspace{10pt}$ $L = a(E(u_2,e) - E(u_1,e))$
   で、得られる。
   B) $0 \lt a \leq b$の場合
   $\hspace{10pt}$ $e = \sqrt{1-(a/b)^2}$
   $\hspace{10pt}$ $ v_1 = th2v(\theta_1)$
   $\hspace{10pt}$ $ v_2 = th2v(\theta_2)$
   とすると、弧長$L$は
   $\hspace{10pt}$ $L = b(E(v_2,e) - E(v_1,e))$
   で、得られる。

5. 計算例
   以下では、角度の単位を数値表の表示に合わせて、「度(°)」とする。
   その2の計算例3)で扱った問題を再度取り挙げる。
   ・シドニー-東京間の子午線距離を求める
   今回も、$E(\varphi,k)$を求めるのに、丸め誤差を防ぐため、専用ソフトを使う。
   ただし、その2で紹介した $ke!^{+}san$ の$E(\phi,k)$は、$\phi$の入力可能範囲が$[-360,360]$となっているが、$[-90,90]$以外では、奇妙な挙動が見られる。$[-90,90]$の範囲だけで利用させて貰うのが、無難と思われる。
   $\hspace{5pt}$ $a = 6378.137km,\,b = 6356.752km$
   $\hspace{5pt}$ $e = \sqrt{1-(6356.752/6378.137)^2} = 0.08181979$
   地理緯度
   　シドニー：南緯33.868333
   　　　東京：北緯35.689556
   地心緯度
   　シドニー：南緯33.690478
   　　　東京：北緯35.507398
   これから、
   $\hspace{5pt}$ $u_1 = th2u(35.507398) = 54.401572$
   $\hspace{5pt}$ $u_2 = th2u(-33.690478) = 123.779349$
   $\hspace{5pt}$ $E(u_2,e) = 2.1559621$
   $\hspace{5pt}$ $E(u_1,e) = 0.9486891$
   よって、
   $\hspace{15pt} L = a(E(u_2,e)-E(u_1,e)) = 7700.15(km)$