fortran：计算第一类椭圆积分

此fortran代码实现了第一类椭圆积分。在此特别向徐士良老先生致敬。本代吗借鉴于徐士良老先生的《fortran常用算法程序集 第二版》。实现代码如下


!//  第一类椭圆积分示例
!// 取k = 0.6, fai = pi/18 * i, i = 0, 1, ..., 10
program main
        implicit none
        integer :: i 
        real*8 :: k, f, melp1, y 
        real*8, parameter :: pi = acos(-1.d0)
        
        k = 0.6
        print*
        write( *, 10 )
10      format( 9x, 'f', 10x, 'FK(k,f)' )

        do i = 0, 10
                f = 0.d0 + real(i,8) * pi/18.d0
                y = melp1(k, f, pi)
                write( *, 30 ) f, y 
        end do
30      format ( 1x, f10.2, 5x, d13.6 )
        print* 


end program main


function melp1(k, f, pi)
        implicit none 
        integer :: n
        real*8 :: k, f, melp1, pi, y, ff, fff, e, fk1

        if ( (k .lt. 0.d0) .or. (k .gt. 1.d0) ) then 
                melp1 = -1.d0+37
                return
        end if

        y = abs(f)
        n = 0

        do
                if ( y .ge. pi ) then
                        n = n + 1
                        y = y - pi 
                else
                        exit
                end  if
        end do

        e = 1.d0 
        if ( y .ge. pi/2.d0 ) then 
                n = n + 1
                e = -1.d0
                y = pi - y 
        end if

        if ( n .eq. 0 ) then 
                ff = fk1( k, y )
                melp1 = ff 
                return 
        else
                ff = fk1(k, pi/2.d0)
                fff = fk1(k, y)
                melp1 = 2.0 * n * ff + e * fff 
                return 
        end  if

end  function melp1 

function fk1(k, y)
        implicit none
        integer :: i, j, m
        real*8 :: t(5), c(5)
        double precision fk1, k, y, a, b, g, s, p, h, aa, bb, w, x, q

        t = [-0.9061798459d0, -0.538469310d0, 0.d0, 0.538469310d0, 0.9061798459d0]
        c = [0.2369268851d0, 0.4786286705d0, 0.5688888889d0, 0.4786286705d0, 0.2369268851d0]

        a = 0.d0
        b = y
        m = 1
        s = (b-a) * 0.02
        p = 0.d0 
        
        do
                h = (b-a) / m 
                g = 0.d0

                do i = 1, m 
                        aa = a + (i-1) * h 
                        bb = a + i * h 
                        w = 0.d0 

                        do j = 1, 5
                                x = ( (bb-aa) * t(j) + (bb+aa) ) / 2.d0
                                fk1 = 1.d0 / sqrt( 1.d0 - k * k * sin(x) * sin(x) )
                                w = w + fk1 * c(j)
                        end do
                
                        g = g + w 

                end do

                g = g * h / 2.d0
                q = abs(g-p) / (1.d0 + abs(g))
        
                if ( (q .ge. 1.d-10) .and. (abs(h) .gt. abs(s)) ) then 
                        p = g 
                        m = m + 1
                else
                        exit
                end if
        
        enddo

        fk1 = g 
        return 
end  function fk1