!=======================================================================
! ch07_adaptive_simpson.f90
!
! Adaptive numerical integration using Simpson's rule.
!
! The routine repeatedly compares Simpson estimates on one interval
! and on its two equal subintervals.  Intervals that do not satisfy
! the local error test are subdivided further.
!
! The algorithm is non-recursive: interval data are stored explicitly
! in arrays that act as a stack.
!
! Example:
!   Integral of
!
!             1
!     -----------------
!     1 + 400(x-0.2)^2
!
!   from 0 to 1.
!
! The integrand has a narrow peak near x = 0.2, making it a useful
! example for adaptive integration.
!
! Exact value:
!
!   [atan(16) + atan(4)] / 20
!
! Original adaptive Simpson routine written by Alexander Godunov,
! July 2012.
! Revised for the companion website, 2026.
!=======================================================================

program main
  implicit none

  double precision :: f
  double precision :: a, b, eps, result, errest, exact
  integer :: nfun
  external :: f

  a   = 0.0d0
  b   = 1.0d0
  eps = 1.0d-8

  exact = (atan(16.0d0) + atan(4.0d0))/20.0d0

  call adaptive_simpson(f, a, b, eps, result, errest, nfun)

  write(*,'(a)')          'Adaptive Simpson integration'
  write(*,'(a)')          'Integral of 1/[1+400(x-0.2)^2] from 0 to 1'
  write(*,'(a,es15.7)')   'Exact value          = ', exact
  write(*,'(a,es15.7)')   'Calculated value     = ', result
  write(*,'(a,es12.4)')   'Estimated error      = ', errest
  write(*,'(a,i10)')      'Function calls       = ', nfun

end program main


double precision function f(x)
  !---------------------------------------------------------------------
  ! Test function for the integration example.
  !
  ! Input:
  !   x - point at which the function is evaluated
  !
  ! Output:
  !   f - value of 1/[1+400(x-0.2)^2]
  !
  ! The narrow peak near x = 0.2 illustrates why adaptive subdivision
  ! is useful: more intervals are introduced where the function varies
  ! rapidly and fewer where it is smooth.
  !---------------------------------------------------------------------
  implicit none

  double precision, intent(in) :: x

  f = 1.0d0/(1.0d0 + 400.0d0*(x-0.2d0)**2)

end function f


subroutine adaptive_simpson(f, a, b, eps, result, errest, nfun)
  !---------------------------------------------------------------------
  ! Adaptive non-recursive integration using Simpson's rule.
  !
  ! Input:
  !   f   - function to integrate
  !   a   - lower integration limit
  !   b   - upper integration limit
  !   eps - requested error tolerance
  !
  ! Output:
  !   result - numerical approximation to the integral
  !   errest - accumulated estimate of the absolute integration error
  !   nfun   - total number of function evaluations
  !
  ! Method:
  !   On each current interval, Simpson's rule is first available for
  !   the whole interval. Two additional function values are then used
  !   to obtain Simpson estimates s1 and s2 on the left and right
  !   halves. The difference
  !
  !       |s1 + s2 - s0|
  !
  !   is used to estimate the local error. For Simpson's rule, the
  !   estimated error in the refined result is approximately this
  !   difference divided by 15.
  !
  !   If the local error test is not satisfied, the interval is divided
  !   into two equal subintervals. Their data are placed on an explicit
  !   stack, so recursion is not required. Function values already
  !   computed at shared points are reused.
  !
  ! Parameter:
  !   im - maximum allowed subdivision depth
  !
  ! Original routine written by Alexander Godunov, July 2012.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: a, b, eps
  double precision, intent(out) :: result, errest
  integer, intent(out) :: nfun
  external :: f

  integer, parameter :: im = 32

  double precision :: tol(im), x(im), h(im)
  double precision :: fa(im), fm(im), fb(im), s(im)
  double precision :: step, err
  double precision :: x0, f0, f1, f2, f3, f4
  double precision :: s0, s1, s2
  integer :: level(im)
  integer :: i, deep

  !---------------------------------------------------------------------
  ! Stage 1: initialize the full interval.
  !
  ! h stores half the width of the current interval. Therefore the
  ! initial Simpson estimate on [a,b] is
  !
  !   h [f(a) + 4 f((a+b)/2) + f(b)] / 3 .
  !
  ! The factor 15 in tol follows from the standard Simpson error
  ! estimate based on the difference between one-panel and two-panel
  ! approximations.
  !---------------------------------------------------------------------

  result = 0.0d0
  errest = 0.0d0
  i = 1

  x(1)     = a
  h(1)     = (b-a)/2.0d0
  fa(1)    = f(a)
  fm(1)    = f(a+h(1))
  fb(1)    = f(b)
  tol(1)   = 15.0d0*eps
  level(1) = 1

  s(1) = h(1)*(fa(1) + 4.0d0*fm(1) + fb(1))/3.0d0
  nfun = 3

  !---------------------------------------------------------------------
  ! Stage 2: process intervals until the stack is empty.
  !---------------------------------------------------------------------

  do while (i > 0)

    ! Evaluate the two new quarter points.
    f1 = f(x(i) +       h(i)/2.0d0)
    f3 = f(x(i) + 3.0d0*h(i)/2.0d0)
    nfun = nfun + 2

    ! Simpson estimates on the left and right halves.
    s1 = h(i)*(fa(i) + 4.0d0*f1 + fm(i))/6.0d0
    s2 = h(i)*(fm(i) + 4.0d0*f3 + fb(i))/6.0d0

    ! Save the current interval before removing it from the stack.
    x0   = x(i)
    f0   = fa(i)
    f2   = fm(i)
    f4   = fb(i)
    step = h(i)
    err  = tol(i)
    s0   = s(i)
    deep = level(i)

    i = i - 1

    ! Accept the refined result when the local error test is satisfied.
    if (abs(s1+s2-s0) <= err) then

      result = result + s1 + s2
      errest = errest + abs(s1+s2-s0)/15.0d0

    else

      if (deep >= im) then
        write(*,'(a)') 'adaptive_simpson: maximum subdivision depth reached.'
        stop
      end if

      ! Right subinterval.
      i = i + 1
      x(i)     = x0 + step
      fa(i)    = f2
      fm(i)    = f3
      fb(i)    = f4
      h(i)     = step/2.0d0
      tol(i)   = err/2.0d0
      s(i)     = s2
      level(i) = deep + 1

      ! Left subinterval.
      i = i + 1
      x(i)     = x0
      fa(i)    = f0
      fm(i)    = f1
      fb(i)    = f2
      h(i)     = h(i-1)
      tol(i)   = tol(i-1)
      s(i)     = s1
      level(i) = level(i-1)

    end if

  end do

end subroutine adaptive_simpson