module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defpi

program f90_CEIS
! Complete Ellipse Integral of the first kind and the second kind
! (Series expansion)
    use defpi
    implicit none
    real(8)::p,Kp,Ep
    integer::ita
    integer::i
    real(4)::t1,t2,ttmax

    t1=DIFT(ttmax)
    do i=0,10
        p=dble(i)*0.1D0
        if(i==10)p=0.9999D0
        call CALKPEP(p,Kp,Ep,ita)
        write(6,'(i5,3(e15.7))') ita,p,Kp,Ep
    end do
    t2=DIFT(ttmax)
    write(6,*) 'time=',t2-t1,'(sec)'

stop

contains

    subroutine CALKPEP(p,Kp,Ep,ita)
        real(8),intent(in)::p
        real(8),intent(out)::Kp,Ep
        integer,intent(out)::ita
        integer::i,j
        real(8)::dK,dE,pn
        real(8),parameter::eps=1.0D-10

        Kp=1.0D0
        Ep=1.0D0
        i=0
        do
            i=i+1
            pn=1.0D0
            do j=1,i
                pn=pn*p*p
            end do
            dK=COEF2(i)*pn
            dE=COEF2(i)*pn/dble(2*i-1)
            Kp=Kp+dK
            Ep=Ep-dE
            if(abs(dK)<eps.and.abs(dE)<eps)exit
            if(100000<i)exit
        end do
        Kp=0.5D0*pi*Kp
        Ep=0.5D0*pi*Ep
        ita=i
    end subroutine CALKPEP

    real(8) function COEF2(n)
        integer,intent(in)::n
        integer::i
        real(8)::f,c1,c2

        f=1.0D0
        do i=1,n
            c1=dble(2*i-1)
            c2=dble(2*i)
            f=f*c1/c2
        end do
        COEF2=f*f
    end function COEF2

    real(4) function DIFT(ttmax)
        real(4),intent(out)::ttmax
        integer::count,count_rate,count_max
        call system_clock(count,count_rate,count_max)
        ttmax=real(count_max)/real(count_rate)
        DIFT=real(count)/real(count_rate)
    end function DIFT

end program f90_CEIS