此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
fortran:计算第一类椭圆积分
最新推荐文章于 2024-08-23 10:38:33 发布