!=======================================================================
! ch07_quanc8.f90
!
! Numerical integration using QUANC8, an adaptive 8-panel
! Newton-Cotes quadrature routine.
!
! Origin:
!   Based on the QUANC8 routine presented by
!   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
!   Computer Methods for Mathematical Computations,
!   Prentice-Hall, 1977.
!
! This version was rewritten and adapted for modern Fortran
! by Alexander Godunov.
!
! Example:
!   Integral of sin(x) from 0 to pi.
!   Exact value = 2.
!
! The program reports the numerical integral, the estimated error,
! the number of function evaluations, and the QUANC8 reliability flag.
! A flag value of zero indicates that the requested tolerance was
! probably satisfied.
!=======================================================================

program ch07_quanc8
  implicit none

  double precision, parameter :: pi = 3.14159265358979323846d0
  double precision :: f, a, b, abserr, relerr
  double precision :: result, errest, flag, exact
  integer :: nofun
  external :: f

  ! Integration interval and requested tolerances
  a      = 0.0d0
  b      = pi
  abserr = 0.0d0
  relerr = 1.0d-8

  call quanc8(f, a, b, abserr, relerr, result, errest, nofun, flag)

  exact = 2.0d0

  write(*,'(a,f16.12)') 'QUANC8 result      = ', result
  write(*,'(a,f16.12)') 'Exact value        = ', exact
  write(*,'(a,es12.4)') 'Actual error       = ', abs(result-exact)
  write(*,'(a,es12.4)') 'Estimated error    = ', errest
  write(*,'(a,i8)')     'Function calls     = ', nofun
  write(*,'(a,f10.4)')  'Reliability flag   = ', flag

end program ch07_quanc8


function f(x)
!-----------------------------------------------------------------------
! Test integrand.
!
! The simple choice f(x) = sin(x) on [0,pi] is used so that the
! numerical result can be checked immediately against the exact value 2.
!-----------------------------------------------------------------------
  implicit none
  double precision :: f, x

  f = sin(x)

end function f


subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag)
!=======================================================================
! Adaptive numerical integration of fun(x) on [a,b].
!
! Method:
!   Automatic adaptive quadrature based on an 8-panel Newton-Cotes rule.
!   The interval is subdivided until the local error estimate satisfies
!   the requested tolerance, subject to limits on refinement and the
!   number of function evaluations.
!
! Input:
!   fun     - integrand function fun(x)
!   a       - lower limit of integration
!   b       - upper limit of integration; b may be less than a
!   abserr  - requested absolute error tolerance (non-negative)
!   relerr  - requested relative error tolerance (non-negative)
!
! Output:
!   result  - numerical approximation to the integral
!   errest  - estimate of the absolute error
!   nofun   - number of integrand evaluations
!   flag    - reliability indicator
!             flag = 0: requested tolerance was probably satisfied
!             flag > 0: one or more convergence limits were encountered
!
! Notes:
!   This is the classic QUANC8 algorithm used here as a self-contained
!   example. The numerical algorithm has been retained; comments and
!   formatting have been cleaned for readability.
!=======================================================================
  implicit none

  double precision :: fun, a, b, abserr, relerr, result, errest, flag
  integer :: nofun

  double precision :: w0, w1, w2, w3, w4
  double precision :: area, x0, f0, stone, step, cor11, temp
  double precision :: qprev, qnow, qdiff, qleft, esterr, tolerr
  double precision :: qright(31), fval(16), x(16)
  double precision :: fsave(8,30), xsave(8,30)
  integer :: levmin, levmax, levout, nomax, nofin
  integer :: lev, nim, i, j
  external :: fun

  ! Stage 1: general initialization and control parameters
  levmin = 1
  levmax = 30
  levout = 6
  nomax  = 5000
  nofin  = nomax - 8*(levmax-levout+2**(levout+1))

  ! Weights for the 8-panel Newton-Cotes rule
  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

  result = 0.0d0
  cor11  = 0.0d0
  errest = 0.0d0
  area   = 0.0d0
  nofun  = 0
  flag   = 0.0d0

  if (a == b) return

  ! Stage 2: initialize the 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
    fval(j) = fun(x(j))
  end do
  nofun = 9

  ! Stage 3: central adaptive calculation
30 continue

  x(1) = (x0+x(2))/2.0d0
  fval(1) = fun(x(1))

  do j = 3, 15, 2
    x(j) = (x(j-1)+x(j+1))/2.0d0
    fval(j) = fun(x(j))
  end do

  nofun = nofun + 8
  step = (x(16)-x0)/16.0d0

  qleft = ( w0*(f0+fval(8))                                      &
          + w1*(fval(1)+fval(7))                                 &
          + w2*(fval(2)+fval(6))                                 &
          + w3*(fval(3)+fval(5))                                 &
          + w4*fval(4) ) * step

  qright(lev+1) = ( w0*(fval(8)+fval(16))                         &
                  + w1*(fval(9)+fval(15))                        &
                  + w2*(fval(10)+fval(14))                       &
                  + w3*(fval(11)+fval(13))                       &
                  + w4*fval(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) go to 50
  if (lev >= levmax) go to 62
  if (nofun > nofin) go to 60
  if (esterr <= tolerr) go to 70

  ! Stage 5: interval has not converged; subdivide it
50 continue
  nim = 2*nim
  lev = lev + 1

  ! Save right-half data for later use
  do i = 1, 8
    fsave(i,lev) = fval(i+8)
    xsave(i,lev) = x(i+8)
  end do

  ! Assemble left-half data for immediate use
  qprev = qleft
  do i = 1, 8
    j = -i
    fval(2*j+18) = fval(j+9)
    x(2*j+18) = x(j+9)
  end do
  go to 30

  ! Stage 6: function-evaluation limit is being approached
60 continue
  nofin  = 2*nofin
  levmax = levout
  flag   = flag + (b-x0)/(b-a)
  go to 70

  ! Maximum refinement level reached
62 continue
  flag = flag + 1.0d0

  ! Stage 7: interval accepted; add its contribution
70 continue
  result = result + qnow
  errest = errest + esterr
  cor11  = cor11 + qdiff/1023.0d0

  ! Locate the next interval
72 continue
  if (nim /= 2*(nim/2)) then
    nim = nim/2
    lev = lev-1
    go to 72
  end if

  nim = nim + 1
  if (lev <= 0) go to 80

  qprev = qright(lev)
  x0    = x(16)
  f0    = fval(16)

  do i = 1, 8
    fval(2*i) = fsave(i,lev)
    x(2*i) = xsave(i,lev)
  end do
  go to 30

  ! Stage 8: finalize the result
80 continue
  result = result + cor11

  ! Ensure that the reported error estimate is not below roundoff level
  if (errest == 0.0d0) return

82 continue
  temp = abs(result) + errest
  if (temp /= abs(result)) return
  errest = 2.0d0*errest
  go to 82

end subroutine quanc8