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

はじめに
　その1では、第I象限に限定して、第2種不完全楕円積分を用いて楕円弧長を求める方法を解説した。今回は、弧の存在範囲を、第I及び第IV象限（すなわち、楕円の右半分）に拡げる。楕円は、上下左右対称なので、第I象限での計算方法が分かれば、基本的には全てのケースで計算可能である。しかし、実際には拡張が必要になった時、短時間で対応する事が困難な事が多く、定式化しておく事は、無意味ではないと思われる。さらに、今回は幾つかの計算例を追加した。

$0 \lt b \leq a$の時の弧長
離心率は、$e = \sqrt{1-(b/a)^2}$である。
その1と同様に、媒介変数$u$を用いた第2の形式で楕円を表すが、楕円の右半分を動くものとする。
$ \hspace{10pt} $ $x = a \sin u,\,y = b \cos u,\,0 \leq u \leq \pi $
$\theta$の動く範囲は、$-\pi/2 \leq \theta \leq \pi/2$とする。
第2種不完全楕円積分$E(u, e)$の計算にあたっては、$\theta$から$u$への変換が必要になるが、分母の0等への注釈の煩わしさを避ける為、変換関数を導入する。
$ \hspace{20pt} th2u(\theta) = \begin{eqnarray} \left\{ \begin{array}{l} \tan^{-1}((b/a)/\tan \theta) &(-\pi/2 \lt \theta \lt \pi/2,\,\theta \neq 0) \\ \pi/2 &(\theta = 0) \\ 0 &(\theta = \pm\pi/2) \end{array} \right. \end{eqnarray} $
また、その1で言及しなかったが、第2種完全楕円積分$E(e)$は四分楕円の弧長に対応している事を指摘しておく。
以下、場合分けして、$[\theta_1,\,\theta_2]$の範囲の弧長$L$を求める。
i) $0 \leq \theta_1 \leq \theta_2 \leq \pi/2$の時
その1と同一条件なので、明らかに、
$\hspace{10pt}$ $L = a(E(u_2,e)-E(u_1,e))$
$\hspace{30pt}$ $u_1 = th2u(\theta_2),\,u_2 = th2u(\theta_1)$
である。
ii) $-\pi/2 \leq \theta_1 \lt 0 \leq \theta_2 \leq \pi/2$の時
弧の範囲を$[\theta_1,\,0]$と$[0,\,\theta_2]$に分けて考察する。
先に$[0,\,\theta_2]$の弧長$L_2$を求める。
$aE(r,e)$は、点$(0,b)$から$u=r$で与えられる点までの弧長であった。$u_1 = th2u(\theta_2)$と置くと、$L_2$は点$(a,0)$から$u=u_1$で与えられる点までの弧長なので、四分楕円の弧長から$aE(u_1,e)$を差し引いて求められる。よって、
$\hspace{10pt}$ $L_2 = a(E(e)-E(u_1,e))$
となる。
次に$[\theta_1,\,0]$の弧長$L_1$だが、楕円は上下対称なので、$u_2 = th2u(-\theta_1)$と置くと、$L_2$と全く同様に扱える。よって、
$\hspace{10pt}$ $L_1 = a(E(e)-E(u_2,e))$
となる。以上から、
$\hspace{10pt}$ $L = L_1 + L_2 = a(2E(e)-E(u_2,e)-E(u_1,e))$
$\hspace{30pt}$ $u_1 = th2u(\theta_2),\,u_2 = th2u(-\theta_1)$
が得られる。
iii) $-\pi/2 \leq \theta_1 \leq \theta_2 \lt 0$の時
楕円の上下対称性を考慮すると明らかに
$\hspace{10pt}$ $L = a(E(u_2,e)-E(u_1,e))$
$\hspace{30pt}$ $u_1 = th2u(-\theta_1),\,u_2 = th2u(-\theta_2)$
となる。

$0 \lt a \leq b$の時の弧長
離心率は、$e = \sqrt{1-(a/b)^2}$である。
その1と同様に、媒介変数$v$を用いた第1の形式で楕円を表すが、楕円の右半分を動くものとする。
$ \hspace{10pt} $ $x = a \cos v,\,y = b \sin v,\,-\pi/2 \leq v \leq \pi/2 $
$\theta$の動く範囲は、2.と同様に$-\pi/2 \leq \theta \leq \pi/2$とする。
また、2.と同様に変換関数を導入する。
$ \hspace{20pt} th2v(\theta) = \begin{eqnarray} \left\{ \begin{array}{l} \tan^{-1}((a/b)\tan \theta) &(-\pi/2 \lt \theta \lt \pi/2) \\ \pm\pi/2 &(\theta = \pm\pi/2) \end{array} \right. \end{eqnarray} $
以下、場合分けして、$[\theta_1,\,\theta_2]$の範囲の弧長$L$を求める。
i) $0 \leq \theta_1 \leq \theta_2 \leq \pi/2$の時
その1と同一条件なので、明らかに、
$\hspace{10pt}$ $L = b(E(v_2,e)-E(v_1,e))$
$\hspace{30pt}$ $v_1 = th2v(\theta_1),\,v_2 = th2v(\theta_2)$
である。
ii) $-\pi/2 \leq \theta_1 \lt 0 \leq \theta_2 \leq \pi/2$の時
弧の範囲を$[\theta_1,\,0]$と$[0,\,\theta_2]$に分けて考察する。
先に$[0,\,\theta_2]$の弧長$L_2$を求める。
$bE(r,e)$は、点$(a,0)$から$v=r$で与えられる点までの弧長であったから$v_2 = th2v(\theta_2)$と置くと、直ちに、
$\hspace{10pt}$ $L_2 = bE(v_2,e)$
となる。
次に$[\theta_1,\,0]$の弧長$L_1$だが、楕円は上下対称なので、$v_1 = th2v(-\theta_1)$と置くと、$L_2$と全く同様に扱える。よって、
$\hspace{10pt}$ $L_1 = bE(v_1,e)$
となる。以上から、
$\hspace{10pt}$ $L = L_1 + L_2 = b(E(v_1,e)+E(v_2,e))$
$\hspace{30pt}$ $v_1 = th2v(-\theta_1),\,v_2 = th2v(\theta_2)$
が得られる。
iii) $-\pi/2 \leq \theta_1 \leq \theta_2 \lt 0$の時
楕円の上下対称性を考慮すると明らかに
$\hspace{10pt}$ $L = b(E(v_2,e)-E(v_1,e))$
$\hspace{30pt}$ $v_1 = th2v(-\theta_2),\,v_2 = th2v(-\theta_1)$
となる。

まとめ
楕円を
$\hspace{10pt}$ $\frac{x^2}{a^2}+ \frac{y^2}{b^2} = 1$
と表す時、$\theta=\theta_1$から$\theta=\theta_2$までの弧長$L$は、第2種不完全楕円積分$E(\varphi, k)$と第2種完全楕円積分$E(k)$を用いて、次のように求められる。ただし、$\theta$は動径と$x$軸との間の反時計方向の角度であり
$\hspace{30pt}$ $-\pi/2 \leq \theta_1 \leq \theta_2 \leq \pi/2$
とする。
A. $0 \lt b \leq a$の場合
$\hspace{5pt}$ $e = \sqrt{1-(b/a)^2}$
$\hspace{5pt}$ i) $0 \leq \theta_1 \leq \theta_2 \leq \pi/2$
$\hspace{10pt}$ $L = a\,(E(u_2,e)-E(u_1,e))$
$\hspace{20pt}$ $u_1 = th2u(\theta_2),\,u_2 = th2u(\theta_1)$
$\hspace{5pt}$ ii) $-\pi/2 \leq \theta_1 \lt 0 \leq \theta_2 \leq \pi/2$
$\hspace{10pt}$ $L = a\,(2E(e)-E(u_2,e)-E(u_1,e))$
$\hspace{20pt}$ $u_1 = th2u(\theta_2),\,u_2 = th2u(-\theta_1)$
$\hspace{5pt}$ iii) $-\pi/2 \leq \theta_1 \leq \theta_2 \lt 0$
$\hspace{10pt}$ $L = a\,(E(u_2,e)-E(u_1,e))$
$\hspace{20pt}$ $u_1 = th2u(-\theta_1),\,u_2 = th2u(-\theta_2)$
B. $0 \lt a \leq b$の場合
$\hspace{5pt}$ $e = \sqrt{1-(a/b)^2}$
$\hspace{5pt}$ i) $0 \leq \theta_1 \leq \theta_2 \leq \pi/2$
$\hspace{10pt}$ $L = b\,(E(v_2,e)-E(v_1,e))$
$\hspace{20pt}$ $v_1 = th2v(\theta_1),\,v_2 = th2v(\theta_2)$
$\hspace{5pt}$ ii) $-\pi/2 \leq \theta_1 \lt 0 \leq \theta_2 \leq \pi/2$
$\hspace{10pt}$ $L = b\,(E(v_2,e)+E(v_1,e))$
$\hspace{20pt}$ $v_1 = th2v(-\theta_1),\,v_2 = th2v(\theta_2)$
$\hspace{5pt}$ iii) $-\pi/2 \leq \theta_1 \leq \theta_2 \lt 0$
$\hspace{10pt}$ $L = b\,(E(v_2,e)-E(v_1,e))$
$\hspace{20pt}$ $v_1 = th2v(-\theta_2),\,v_2 = th2v(-\theta_1)$

計算例
以下では、角度の単位を数値表の表示に合わせて、「度(°)」とする。
1)$a=10,\,b=6$の楕円の、$\theta$=45°の点から点$(0,b)$までの弧長を求める
ケースA)-i)に該当する。
$e = \sqrt{1-(6/10)^2} = 0.8$
$u_1$は計算するまでもなく0
$u_2 = th2u(45) = 30.96376$
$E(31,0.8)= 0.52468$
よって、
$\hspace{20pt} L = 10E(31,0.8) = 5.2468$
2)$a=100,\,b=125$の楕円の、$\theta$=15°の点から$\theta$=60°の点までの弧長を求める
ケースB)-i)に該当する。
$e = \sqrt{1-(100/125)^2} = 0.6$
$v_1 = th2v(15) = 12.09879$
$v_2 = th2v(60) = 54.18247$
$E(12,0.6)= 0.20889$
$E(54,0.6)= 0.89872$
よって、
$\hspace{20pt} L = 125(E(54,0.6)-E(12,0.6)) = 86.2$
3)シドニー-東京間の子午線距離を求める
（任意の子午線上の、シドニーの緯度の点から、東京の緯度の点までの距離を求める例で、シドニー-東京間の距離を求めるのではない。念の為。）
ここでは、$E(\varphi,k)$を求めるのに、丸め誤差を防ぐため、専用ソフトを使う。
ネット上では、たとえば $ke!^{+}san$ が利用できる。
a = 6378.137km, b = 6356.752km
シドニー：南緯33.868333
　　東京：北緯35.689556
ケースA)-ii)に該当する。
$e = \sqrt{1-(6356.752/6378.137)^2} = 0.08181979$
ここで、$u_1,u_2$を求めるのだが、上記緯度は、詳しくは地理緯度と呼ばれるもので、楕円上の点に立てた法線が赤道面となす角度になっている。（詳細は 緯度 を参照）$\theta$として使えるのは地球中心を通る地心緯度なので、先に地理緯度から地心緯度への変換が必要である。地理緯度を$\theta'$、地心緯度を$\theta$とすると、変換式は次のようになる。
$\theta = \tan^{-1}((1-e^2)\tan \theta')$
これから、地心緯度は
シドニー：南緯33.690478
　　東京：北緯35.507398
となる。
$u_1 = th2u(35.507398) = 54.401572$
$u_2 = th2u(33.690478) = 56.220652$
$E(u_2,e) = 0.9803661$
$E(u_1,e) = 0.9486891$
$E(e) = 1.5681641$
よって、
$\hspace{10pt} L = a(2E(e)-E(u_2,e)-E(u_1,e)) = 7700.15(km)$
この値は国土地理院の測地線に沿って2点間の距離を出してくれる 測量計算 サイトでの結果と完全に一致している。