!=======================================================================
! ch07_adaptive_gauss_kronrod.f90
!
! Adaptive numerical integration using embedded Gauss-Kronrod rules.
!
! Two quadrature pairs are included:
!   1. Gauss 7-point / Kronrod 15-point (G7-K15)
!   2. Gauss 10-point / Kronrod 21-point (G10-K21)
!
! The same adaptive subdivision algorithm is used with either pair.
!
! Example:
!   Integral of 1/[1+400(x-0.2)^2] from 0 to 1.
!
! Exact value:
!   [atan(16) + atan(4)] / 20
!
! The integrand has a narrow peak near x = 0.2, making it a useful
! example for adaptive integration.
!
! Written by Alexander Godunov.
! Revised for the companion website, 2026.
!=======================================================================

program main
  implicit none

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

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

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

  write(*,'(a)') 'Adaptive Gauss-Kronrod integration'
  write(*,'(a)') 'Integral of 1/[1+400(x-0.2)^2] from 0 to 1'
  write(*,'(a,es15.7)') 'Exact value          = ', exact
  write(*,*)

  call GKinteg(f, a, b, eps, result, errest, nfun, 1)
  write(*,'(a,es15.7)') 'G7-K15 result        = ', result
  write(*,'(a,es12.4)') 'Estimated error      = ', errest
  write(*,'(a,i10)')    'Function calls       = ', nfun
  write(*,*)

  call GKinteg(f, a, b, eps, result, errest, nfun, 2)
  write(*,'(a,es15.7)') 'G10-K21 result       = ', result
  write(*,'(a,es12.4)') 'Estimated error      = ', errest
  write(*,'(a,i10)')    'Function calls       = ', nfun

end program main


double precision function f(x)
  !---------------------------------------------------------------------
  ! Test integrand: a narrow peak centered near x = 0.2.
  !---------------------------------------------------------------------
  implicit none
  double precision, intent(in) :: x

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

end function f


subroutine GKinteg(fun, a, b, eps, result, errest, nfun, key)
  !---------------------------------------------------------------------
  ! Adaptive non-recursive Gauss-Kronrod integration.
  !
  ! Input:
  !   fun  - function to integrate
  !   a,b  - lower and upper limits
  !   eps  - requested tolerance used in the local convergence test
  !   key  - quadrature pair:
  !            1 : Gauss 7  / Kronrod 15
  !            2 : Gauss 10 / Kronrod 21
  !
  ! Output:
  !   result - numerical value of the integral
  !   errest - accumulated Gauss-Kronrod difference
  !   nfun   - number of function evaluations
  !---------------------------------------------------------------------
  implicit none

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

  integer, parameter :: im = 32
  double precision :: tol(im), x(im), h(im)
  double precision :: x0, step, err
  double precision :: gauss, kronrod
  integer :: level(im)
  integer :: i, deep
  external :: fun

  result = 0.0d0
  errest = 0.0d0
  nfun   = 0
  i      = 1

  x(1)     = a
  h(1)     = b-a
  tol(1)   = eps
  level(1) = 1

  do while (i > 0)

    select case (key)
      case (1)
        call G7K15(fun, x(i), x(i)+h(i), gauss, kronrod)
        nfun = nfun + 15
      case (2)
        call G10K21(fun, x(i), x(i)+h(i), gauss, kronrod)
        nfun = nfun + 21
      case default
        write(*,'(a)') 'GKinteg error: key must be 1 or 2.'
        stop
    end select

    x0   = x(i)
    step = h(i)
    err  = tol(i)
    deep = level(i)

    i = i - 1

    if ((200.0d0*abs(kronrod-gauss))**1.5d0 <= err) then
      result = result + kronrod
      errest = errest + abs(kronrod-gauss)
    else
      if (deep >= im) then
        write(*,'(a)') 'GKinteg error: maximum subdivision depth reached.'
        stop
      end if

      ! Right subinterval
      i = i + 1
      h(i)     = step/2.0d0
      tol(i)   = err/2.0d0
      x(i)     = x0 + h(i)
      level(i) = deep + 1

      ! Left subinterval
      i = i + 1
      x(i)     = x0
      h(i)     = h(i-1)
      tol(i)   = tol(i-1)
      level(i) = level(i-1)
    end if

  end do

end subroutine GKinteg

subroutine G7K15(fun, a, b, g7, k15)
  !---------------------------------------------------------------------
  ! Apply the Gauss 7-point and Kronrod 15-point rules on [a,b].
  ! The Gauss points are a subset of the Kronrod points.
  !---------------------------------------------------------------------
  implicit none

  double precision :: fun
  double precision, intent(in)  :: a, b
  double precision, intent(out) :: g7, k15
  double precision :: h, c, f0
  double precision :: x(0:7), wg(0:3), wk(0:7)
  double precision :: fl(7), fr(7)
  integer :: i
  external :: fun

  data x /0.000000000000000d0, 0.207784955007898d0, &
          0.405845151377397d0, 0.586087235467691d0, &
          0.741531185599394d0, 0.864864423359769d0, &
          0.949107912342758d0, 0.991455371120813d0/

  data wg /0.417959183673469d0, 0.381830050505119d0, &
           0.279705391489277d0, 0.129484966168870d0/

  data wk /0.209482141084728d0, 0.204432940075298d0, &
           0.190350578064785d0, 0.169004726639267d0, &
           0.140653259715525d0, 0.104790010322250d0, &
           0.063092092629979d0, 0.022935322010529d0/

  h = (b-a)/2.0d0
  c = (b+a)/2.0d0

  f0 = fun(c)
  do i = 1, 7
    fl(i) = fun(c - h*x(i))
    fr(i) = fun(c + h*x(i))
  end do

  g7 = wg(0)*f0
  do i = 1, 3
    g7 = g7 + wg(i)*(fl(2*i) + fr(2*i))
  end do
  g7 = h*g7

  k15 = wk(0)*f0
  do i = 1, 7
    k15 = k15 + wk(i)*(fl(i) + fr(i))
  end do
  k15 = h*k15

end subroutine G7K15


subroutine G10K21(fun, a, b, g10, k21)
  !---------------------------------------------------------------------
  ! Apply the Gauss 10-point and Kronrod 21-point rules on [a,b].
  !---------------------------------------------------------------------
  implicit none

  double precision :: fun
  double precision, intent(in)  :: a, b
  double precision, intent(out) :: g10, k21
  double precision :: h, c, f0
  double precision :: x(11), wg(5), wk(11)
  double precision :: fl(11), fr(11)
  integer :: i
  external :: fun

  data x /0.0000000000000000000000000d0, &
          0.1488743389816312108848260d0, &
          0.2943928627014601981311266d0, &
          0.4333953941292471907992659d0, &
          0.5627571346686046833390001d0, &
          0.6794095682990244062343274d0, &
          0.7808177265864168970637176d0, &
          0.8650633666889845107320967d0, &
          0.9301574913557082260012072d0, &
          0.9739065285171717200779640d0, &
          0.9956571630258080807355273d0/

  data wk /0.1494455540029169056649365d0, &
           0.1477391049013384913748415d0, &
           0.1427759385770600807970943d0, &
           0.1347092173114733259280540d0, &
           0.1234919762620658510779581d0, &
           0.1093871588022976418992106d0, &
           0.0931254545836976055350655d0, &
           0.0750396748109199527670431d0, &
           0.0547558965743519960313813d0, &
           0.0325581623079647274788190d0, &
           0.0116946388673718742780644d0/

  data wg /0.2955242247147528701738930d0, &
           0.2692667193099963550912269d0, &
           0.2190863625159820439955349d0, &
           0.1494513491505805931457763d0, &
           0.0666713443086881375935688d0/

  h = (b-a)/2.0d0
  c = (b+a)/2.0d0

  f0 = fun(c)
  do i = 2, 11
    fl(i) = fun(c - h*x(i))
    fr(i) = fun(c + h*x(i))
  end do

  g10 = 0.0d0
  do i = 1, 5
    g10 = g10 + wg(i)*(fl(2*i) + fr(2*i))
  end do
  g10 = h*g10

  k21 = wk(1)*f0
  do i = 2, 11
    k21 = k21 + wk(i)*(fl(i) + fr(i))
  end do
  k21 = h*k21

end subroutine G10K21