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

program f90_CEIG
! 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 CEI_Gauss(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

    real(8) function CALFKP(x,p)
        real(8),intent(in)::x,p
        CALFKP=1.0D0/sqrt(1.0D0-p*p*sin(x)*sin(x))
    end function CALFKP

    real(8) function CALFEP(x,p)
        real(8),intent(in)::x,p
        CALFEP=sqrt(1.0D0-p*p*sin(x)*sin(x))
    end function CALFEP

    subroutine CEI_Gauss(p,sKP,sEP,ita)
        real(8),intent(in)::p
        real(8),intent(out)::sKP,sEP
        integer,intent(out)::ita
        integer,parameter::nt=10
        integer,parameter::n=10
        integer,parameter::mmax=1000
        real(8),parameter::a=0.0D0
        real(8),parameter::b=0.5D0*pi
        real(8),parameter::eps=1.0D-10
        integer::i,j,k
        real(8),allocatable::gu(:),gw(:)
        real(8)::x,fKP,fEP
        real(8)::sKPi,sEPi,ai,bi,sKP0,sEP0

        allocate(gu(1:nt))
        allocate(gw(1:nt))

        call GUW(n,gu,gw,nt)
        sKP0=0.0D0
        sEP0=0.0D0
        k=0
        do
            k=k+1
            sKP=0.0D0
            sEP=0.0D0
            do j=1,k
                ai=a+(b-a)/dble(k)*dble(j-1)
                bi=a+(b-a)/dble(k)*dble(j)
                sKPi=0.0D0
                sEPi=0.0D0
                do i=1,n
                    x=0.5*(ai+bi)+0.5*(bi-ai)*gu(i)
                    fKP=CALFKP(x,p)
                    fEP=CALFEP(x,p)
                    sKPi=sKPi+gw(i)*fKP
                    sEPi=sEPi+gw(i)*fEP
                end do
                sKPi=0.5*(bi-ai)*sKPi
                sEPi=0.5*(bi-ai)*sEPi
                sKP=sKP+sKPi
                sEP=sEP+sEPi
            end do
            if(abs(sKP0-sKP)<eps.and.abs(sEP0-sEP)<eps)exit
            sKP0=sKP
            sEP0=sEP
            if(mmax<k)exit
        end do
        ita=k
    end subroutine CEI_Gauss

    subroutine GUW(n,gu,gw,nt)
        integer,intent(in)::n,nt
        real(8),intent(out)::gu(nt),gw(nt)

        select case(n)
        case(2)
            gu( 1)= -0.577350269189625764509D0; gw( 1)=  1.000000000000000000000D0
            gu( 2)=  0.577350269189625764509D0; gw( 2)=  1.000000000000000000000D0
        case(5)
            gu( 1)= -0.906179845938663992798D0; gw( 1)=  0.236926885056189087514D0
            gu( 2)= -0.538469310105683091036D0; gw( 2)=  0.478628670499366468041D0
            gu( 3)=  0.000000000000000000000D0; gw( 3)=  0.568888888888888888889D0
            gu( 4)=  0.538469310105683091036D0; gw( 4)=  0.478628670499366468041D0
            gu( 5)=  0.906179845938663992798D0; gw( 5)=  0.236926885056189087514D0
        case(10)
            gu( 1)= -0.973906528517171720078D0; gw( 1)=  0.066671344308688137594D0
            gu( 2)= -0.865063366688984510732D0; gw( 2)=  0.149451349150580593146D0
            gu( 3)= -0.679409568299024406234D0; gw( 3)=  0.219086362515982043996D0
            gu( 4)= -0.433395394129247190799D0; gw( 4)=  0.269266719309996355091D0
            gu( 5)= -0.148874338981631210885D0; gw( 5)=  0.295524224714752870174D0
            gu( 6)=  0.148874338981631210885D0; gw( 6)=  0.295524224714752870174D0
            gu( 7)=  0.433395394129247190799D0; gw( 7)=  0.269266719309996355091D0
            gu( 8)=  0.679409568299024406234D0; gw( 8)=  0.219086362515982043996D0
            gu( 9)=  0.865063366688984510732D0; gw( 9)=  0.149451349150580593146D0
            gu(10)=  0.973906528517171720078D0; gw(10)=  0.066671344308688137594D0
        end select
    end subroutine GUW

    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_CEIG