下面代码的功能实现了第二类椭圆积分的计算。
本代码借鉴于徐士良老先生的《fortran常用算法程序集 第二版》,再次向老先生致敬!
!// 第二类椭圆积分示例
!// 取k = 0.5, fai = pi/18 * i, i = 0, 1, ..., 10
program main
implicit none
integer :: i
real*8 :: k, f, y, melp2
real*8, parameter :: pi = acos(-1.d0)
k = 0.5
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 = melp2(k, f, pi)
write( *, 30 ) f, y
end do
30 format ( 1x, f10.2, 5x, d13.6 )
print*
end program main
function melp2(k, f, pi)
implicit none
integer :: n
real*8 :: k, f, melp2, pi, y, ff, fff, e, ek1
if ( (k .lt. 0.d0) .or. (k .gt. 1.d0) ) then
melp2 = -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 = ek1( k, y )
melp2 = ff
return
else
ff = ek1(k, pi/2.d0)
fff = ek1(k, y)
melp2 = 2.0 * n * ff + e * fff
return
end if
end function melp2
function ek1(k, y)
implicit none
integer :: i, j, m
real*8 :: t(5), c(5)
real*8 :: ek1, 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
ek1 = sqrt( 1.d0 - k * k * sin(x) * sin(x) )
w = w + ek1 * 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
ek1 = g
return
end function ek1