subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag)
!! Estimate the integral of fun(x) from a to b
!! author: P J Knight, CCFE, Culham Science Centre
!! fun : external function : integrand function subprogram fun(x)
!! a : input real : lower limit of integration
!! b : input real : upper limit of integration (b may be less than a)
!! abserr : input real : absolute error tolerance (should be non-negative)
!! relerr : input real : relative error tolerance (should be non-negative)
!! result : output real : approximation to the integral hopefully
!! satisfying the least stringent of the two error tolerances
!! errest : output real : estimate of the magnitude of the actual error
!! nofun : output integer : number of function values used in calculation
!! flag : output real : Reliability indicator; if flag is zero, then
!! result probably satisfies the error tolerance. If flag is
!! xxx.yyy , then xxx = the number of intervals which have
!! not converged and 0.yyy = the fraction of the interval
!! left to do when the limit on nofun was approached.
!! This routine estimates the integral of fun(x) from a to b
!! to a user provided tolerance. An automatic adaptive
!! routine based on the 8-panel Newton-Cotes rule.
!! http://www.netlib.org/fmm/index.html :
!! Computer Methods for Mathematical Computations,
!! G E Forsythe, M A Malcolm, and C B Moler,
!! Prentice-Hall, Englewood Cliffs, New Jersey
!! 1977, ISBN 0-13-165332-6
!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
interface
function fun(rho) result(fint)
#ifndef dp
use, intrinsic :: iso_fortran_env, only: dp=>real64
#endif
real(dp), intent(in) :: rho
real(dp) :: fint
end function fun
end interface
! Arguments
external :: fun
real(dp), intent(in) :: a, b, abserr, relerr
real(dp), intent(out) :: result, errest, flag
integer, intent(out) :: nofun
! Local variables
real(dp) :: w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp
real(dp) :: qprev,qnow,qdiff,qleft,esterr,tolerr
real(dp), dimension(31) :: qright
real(dp), dimension(16) :: f, x
real(dp), dimension(8,30) :: fsave, xsave
integer :: levmin,levmax,levout,nomax,nofin,lev,nim,i,j
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Stage 1: general initialization
! Set constants
levmin = 1
levmax = 30
levout = 6
nomax = 5000
nofin = nomax - 8*(levmax-levout+2**(levout+1))
! Trouble when nofun reaches nofin
w0 = 3956.0D0 / 14175.0D0
w1 = 23552.0D0 / 14175.0D0
w2 = -3712.0D0 / 14175.0D0
w3 = 41984.0D0 / 14175.0D0
w4 = -18160.0D0 / 14175.0D0
! Initialize running sums to zero
flag = 0.0D0
result = 0.0D0
cor11 = 0.0D0
errest = 0.0D0
area = 0.0D0
nofun = 0
if (a == b) return
! Stage 2: initialization for first interval
lev = 0
nim = 1
x0 = a
x(16) = b
qprev = 0.0D0
f0 = fun(x0)
stone = (b - a) / 16.0D0
x(8) = (x0 + x(16)) / 2.0D0
x(4) = (x0 + x(8)) / 2.0D0
x(12) = (x(8) + x(16)) / 2.0D0
x(2) = (x0 + x(4)) / 2.0D0
x(6) = (x(4) + x(8)) / 2.0D0
x(10) = (x(8) + x(12)) / 2.0D0
x(14) = (x(12) + x(16)) / 2.0D0
do j = 2, 16, 2
f(j) = fun(x(j))
end do
nofun = 9
! Stage 3: central calculation
! Requires qprev,x0,x2,x4,...,x16,f0,f2,f4,...,f16
! Calculates x1,x3,...x15, f1,f3,...f15,qleft,qright,qnow,qdiff,area
main_loop: do
x(1) = (x0 + x(2)) / 2.0D0
f(1) = fun(x(1))
do j = 3, 15, 2
x(j) = (x(j-1) + x(j+1)) / 2.0D0
f(j) = fun(x(j))
end do
nofun = nofun + 8
step = (x(16) - x0) / 16.0D0
qleft = (w0*(f0 + f(8)) + w1*(f(1)+f(7)) + w2*(f(2)+f(6)) &
+ w3*(f(3)+f(5)) + w4*f(4)) * step
qright(lev+1) = (w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14)) &
+ w3*(f(11)+f(13)) + w4*f(12)) * step
qnow = qleft + qright(lev+1)
qdiff = qnow - qprev
area = area + qdiff
! Stage 4: interval convergence test
esterr = abs(qdiff) / 1023.0D0
tolerr = max(abserr,relerr*abs(area)) * (step/stone)
if ( (lev < levmin).or. &
((lev < levmax).and.(nofun <= nofin) &
.and.(esterr > tolerr)) ) then
! Stage 5: no convergence
! Locate next interval
nim = 2*nim
lev = lev+1
! Store right hand elements for future use
do i = 1, 8
fsave(i,lev) = f(i+8)
xsave(i,lev) = x(i+8)
end do
! Assemble left hand elements for immediate use
qprev = qleft
do i = 1, 8
j = -i
f(2*j+18) = f(j+9)
x(2*j+18) = x(j+9)
end do
cycle main_loop
else if (lev >= levmax) then
flag = flag + 1.0D0
else if (nofun > nofin) then
! Stage 6: trouble section
! Number of function values is about to exceed limit
nofin = 2*nofin
levmax = levout
flag = flag + (b - x0) / (b - a)
end if
! Stage 7: interval converged
! Add contributions into running sums
result = result + qnow
errest = errest + esterr
cor11 = cor11 + qdiff / 1023.0D0
! Locate next interval
do
if (nim == (2*(nim/2))) exit
nim = nim/2
lev = lev-1
end do
nim = nim + 1
if (lev <= 0) exit main_loop
! Assemble elements required for the next interval
qprev = qright(lev)
x0 = x(16)
f0 = f(16)
do i = 1, 8
f(2*i) = fsave(i,lev)
x(2*i) = xsave(i,lev)
end do
end do main_loop
! Stage 8: finalize and return
result = result + cor11
! Make sure errest not less than roundoff level
if (errest == 0.0D0) return
estimate_error: do
temp = abs(result) + errest
if (temp /= abs(result)) exit estimate_error
errest = 2.0D0*errest
end do estimate_error
end subroutine quanc8