0

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

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

  2. $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)$
    となる。

  3. $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)$
    となる。

  4. まとめ
    楕円を
    $\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)$

  5. 計算例
    以下では、角度の単位を数値表の表示に合わせて、「度(°)」とする。
    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点間の距離を出してくれる 測量計算 サイトでの結果と完全に一致している。

投稿日:2022423

この記事を高評価した人

高評価したユーザはいません

この記事に送られたバッジ

バッジはありません。

投稿者

GonJii
2
1258

コメント

他の人のコメント

コメントはありません。
読み込み中...
読み込み中