第一種完全楕円積分 [fortran]

こんにちは　管理人のさくらえびです。

本稿では、第1種完全楕円積分のプログラムを紹介します。
自分で投稿したつもりでアップロードできていなかったのかも・・・
えー、はい。英語で言うと、
complete elliptic Integral of the 1st kind
さて、まずは数式のご紹介
下記定積分を第1種完全楕円積分と言います。
\[
K(k) = \int_0^\frac{\pi}{2} \frac{d\theta}{(1-k^2sin^2\theta)^{1/2}}
\]
プログラムしやすい式に書き直すと、
\[
K(k) = \frac{\pi}{2}\sum_{n=0}^\infty \bigl(\frac{(2n-1)!!}{(2n)!!} \bigr)^2 {k^{2n}}
\]
こんな感じ。
でも、収束するまで時間がかかるので、
\[
\kappa = \frac{1-\sqrt{1-k^2}}{1+\sqrt{1-k^2}}
\]
と、ベキ級数にしてしまえば、
\[
E(k) = \frac{\pi (1+\kappa)}{2}\bigl[1+\sum_{n=1}^\infty \frac{(2n-1)!!}{(2n)!!}\kappa^{2n}\bigr]
\]
こんな感じです。
このプログラミングは下記のとおりです。

function ELLIP_1st(x,endT)
!*----------------------------------------------------
!   This function derives Complete Elliptic Integral of the 1st kind
!-----------------------------------------------------
!definision valiables
implicit none
double precision ,parameter :: pi = 3.141592653589793d0
double precision ELLIP_1st,tempE
double precision x,xx,FACT2,CC,CCE
integer*4 i,endT
xx = (1.d0-dsqrt(1.d0-x**2))/(1.d0+dsqrt(1.d0-x**2))
CCE = pi*0.5d0*(1.d0+xx)
tempE = 1
do i = 1 ,endT
CC = (FACT2(2*i-1)/FACT2(2*i))**2
tempE = tempE + CC*+xx**(2.d0*dble(i))
enddo
ELLIP_1st = CCE*tempE
end function ELLIP_1st

階乗プログラムについてはこちら
ちなみに数表は、下記の本を参考にデバッグしております。
50回転くらいしてみると、こんな結果が出ました。
k = 1で無限になるのですが、飛び方が弱いですね。どうせ1は使わないので、こんなところでしょう。

これで、完全楕円積分は1種と2種揃いました。第二種完全楕円積分はこちら