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

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

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

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

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

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

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

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