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)
        real*8 :: fk1, k, y, a, b, g, s, p, h, aa, bb, w, x, q

        !// 5个高斯零点与权系数
        t = [-0.9061798459d0, -0.538469310d0, 0.d0, 0.538469310d0, 0.9061798459d0]
        c = [0.2369268851d0, 0.4786286705d0, 0.5688888889d0, 0.4786286705d0, 0.2369268851d0]


        // 10个高斯零点与权系数
        !t = [-0.973906528517172d0,-0.865063366688985d0,-0.679409568299024d0,-0.433395394129247d0,-0.148874338981631d0,0.148874338981631d0,0.433395394129247d0,0.679409568299024d0,0.865063366688985d0,0.973906528517172d0]
        !c = [6.667134430868699d-2, 0.149451349150581d0, 0.219086362515982d0, 0.269266719309996d0, 0.295524224714753d0,0.295524224714753d0,0.269266719309996d0,0.219086362515982d0,0.149451349150581d0,6.667134430868939d-2]
 
        !// 20个高斯零点与权系数
        !t = [-0.993128599185095d0,      -0.963971927277914d0,      -0.912234428251326d0,&
        !     -0.839116971822219d0,      -0.746331906460151d0,      -0.636053680726515d0,&
        !     -0.510867001950827d0,      -0.373706088715420d0,      -0.227785851141645d0,&
        !     -7.652652113349734d-2,      7.652652113349734d-2,      0.227785851141645d0,&
        !      0.373706088715420d0,       0.510867001950827d0,       0.636053680726515d0,&
        !      0.746331906460151d0,       0.839116971822219d0,       0.912234428251326d0,&
        !      0.963971927277914d0,       0.993128599185095d0]
        !c = [1.761400713915310d-2,       4.060142980038731d-2,      6.267204833410855d-2,&
        !     8.327674157670428d-2,       0.101930119817240d0,       0.118194531961518d0,&
        !     0.131688638449176d0,        0.142096109318382d0,       0.149172986472604d0,&
        !     0.152753387130726d0,        0.152753387130726d0,       0.149172986472604d0,&
        !     0.142096109318382d0,        0.131688638449176d0,       0.118194531961518d0,&
        !     0.101930119817240d0,        8.327674157670428d-2,      6.267204833410855d-2,&
        !     4.060142980038731d-2,       1.761400713915310d-2]
        
        !// 50个高斯零点与权系数
        !t = [ -0.998866404420071d0,      -0.994031969432091d0,      -0.985354084048006d0,&
        !      -0.972864385106692d0,      -0.956610955242808d0,      -0.936656618944878d0,&
        !      -0.913078556655792d0,      -0.885967979523613d0,      -0.855429769429946d0,&
        !      -0.821582070859336d0,      -0.784555832900399d0,      -0.744494302226068d0,&
        !      -0.701552468706822d0,      -0.655896465685439d0,      -0.607702927184950d0,&
        !      -0.557158304514650d0,      -0.504458144907464d0,      -0.449806334974039d0,&
        !      -0.393414311897565d0,      -0.335500245419437d0,      -0.276288193779532d0,&
        !      -0.216007236876042d0,      -0.154890589998146d0,      -9.317470156008614d-2,&
        !      -3.109833832718887d-2,      3.109833832718888d-2,      9.317470156008614d-2,&
        !       0.154890589998146d0,       0.216007236876042d0,       0.276288193779532d0,&
        !       0.335500245419437d0,       0.393414311897565d0,       0.449806334974039d0,&
        !       0.504458144907464d0,       0.557158304514650d0,       0.607702927184950d0,&
        !       0.655896465685439d0,       0.701552468706822d0,       0.744494302226068d0,&
        !       0.784555832900399d0,       0.821582070859336d0,       0.855429769429946d0,&
        !       0.885967979523613d0,       0.913078556655792d0,       0.936656618944878d0,&
        !       0.956610955242808d0,       0.972864385106692d0,       0.985354084048006d0,&
        !       0.994031969432091d0,       0.998866404420071d0]
        !
        !c = [2.908622553156514d-003,  6.759799195745833d-003,  1.059054838365176d-002,&
        !     1.438082276148574d-002,  1.811556071348932d-002,  2.178024317012469d-002,&
        !     2.536067357001282d-002,  2.884299358053498d-002,  3.221372822357788d-002,&
        !     3.545983561514593d-002,  3.856875661258736d-002,  4.152846309014788d-002,&
        !     4.432750433880345d-002,  4.695505130394848d-002,  4.940093844946638d-002,&
        !     5.165570306958110d-002,  5.371062188899625d-002,  5.555774480621256d-002,&
        !     5.718992564772835d-002,  5.860084981322244d-002,  5.978505870426543d-002,&
        !     6.073797084177023d-002,  6.145589959031672d-002,  6.193606742068323d-002,&
        !     6.217661665534728d-002,  6.217661665534728d-002,  6.193606742068323d-002,&
        !     6.145589959031672d-002,  6.073797084177023d-002,  5.978505870426543d-002,&
        !     5.860084981322244d-002,  5.718992564772835d-002,  5.555774480621256d-002,&
        !     5.371062188899625d-002,  5.165570306958110d-002,  4.940093844946638d-002,&
        !     4.695505130394848d-002,  4.432750433880345d-002,  4.152846309014788d-002,&
        !     3.856875661258789d-002,  3.545983561514593d-002,  3.221372822357788d-002,&
        !     2.884299358053498d-002,  2.536067357001282d-002,  2.178024317012469d-002,&
        !     1.811556071348932d-002,  1.438082276148574d-002,  1.059054838365176d-002,&
        !     6.759799195745833d-003,  2.908622553156514d-003]


        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