!=======================================================================
! ch07_nested_simpson_2d.f90
!
! Two-dimensional integration using nested Simpson integration.
!
! The program evaluates
!
!          1   sin(x)
!          /      /
!         |      |       x^2
!     I = |      |   --------- dy dx .
!         |      |    y^2 + 2
!         /      /
!          0      0
!
! The two-dimensional problem is reduced to repeated one-dimensional
! integrations:
!
!              sin(x)
!                 /
!        F(x) =  |  f(x,y) dy,
!                 /
!                 0
!
!        I = integral F(x) dx,   0 <= x <= 1.
!
! The same Simpson routine is used for both the inner and outer
! integrations.  The number of intervals is doubled until
!
!        |I_(2n) - I_n| / 15 < tolerance.
!
! This is automatic global refinement rather than local adaptive
! subdivision.  It illustrates the nested-integration strategy
! discussed in Section 7.7, "Multidimensional integration."
!
! Based on code written by Alexander Godunov, October 2009.
! Revised and reorganized for the companion website, 2026.
!=======================================================================

module nested_data
  implicit none

  ! Current value of x used by the inner y integration.
  double precision :: x_current

  ! Tolerance used for each inner integral.
  double precision :: eps_inner

  ! Count evaluations of the original two-dimensional integrand.
  integer :: nfun_2d

end module nested_data


program main
  use nested_data, only: eps_inner, nfun_2d
  implicit none

  double precision :: inner_integral
  double precision :: a, b
  double precision :: eps_outer
  double precision :: integral, error_estimate
  double precision :: reference, absolute_error
  integer :: n_used
  external :: inner_integral

  a = 0.0d0
  b = 1.0d0

  ! Tolerances for the nested calculation.
  ! The inner tolerance is chosen somewhat smaller than the outer one
  ! so that errors in F(x) do not dominate the outer integration.
  eps_outer = 1.0d-7
  eps_inner = 1.0d-8

  ! High-accuracy reference value for the test problem.
  ! If the integrand or limits are changed, replace or omit this value.
  reference = 1.0344649764293148d-1

  nfun_2d = 0

  call simpson_refine(inner_integral, a, b, eps_outer,              &
                      integral, error_estimate, n_used)

  absolute_error = abs(integral-reference)

  write(*,'(a)')        'Two-dimensional nested Simpson integration'
  write(*,'(a)')        'Integral of x^2/(y^2+2), 0 <= x <= 1,'
  write(*,'(a)')        'with 0 <= y <= sin(x)'
  write(*,'(a,es15.7)') 'Reference value      = ', reference
  write(*,'(a,es15.7)') 'Calculated value     = ', integral
  write(*,'(a,es12.4)') 'Estimated outer error= ', error_estimate
  write(*,'(a,es12.4)') 'Absolute error       = ', absolute_error
  write(*,'(a,i10)')    'Outer intervals      = ', n_used
  write(*,'(a,i10)')    '2D function calls    = ', nfun_2d

end program main


double precision function integrand(x,y)
  !---------------------------------------------------------------------
  ! Original two-dimensional function f(x,y).
  !
  ! To use this program for another problem, replace only the line
  ! defining integrand(x,y) below.  The nested Simpson routine does not
  ! need to be changed.
  !
  ! Input:
  !   x, y - integration variables
  !
  ! Output:
  !   integrand - value of f(x,y)
  !---------------------------------------------------------------------
  use nested_data, only: nfun_2d
  implicit none

  double precision, intent(in) :: x, y

  nfun_2d = nfun_2d + 1

  integrand = x*x/(y*y + 2.0d0)

end function integrand


double precision function y_lower(x)
  !---------------------------------------------------------------------
  ! Lower limit of the inner y integration.
  !
  ! Change this function if the lower boundary depends on x.
  !---------------------------------------------------------------------
  implicit none

  double precision, intent(in) :: x

  ! Multiplication by x keeps the dependence explicit while the
  ! present lower boundary remains identically zero.
  y_lower = 0.0d0*x

end function y_lower


double precision function y_upper(x)
  !---------------------------------------------------------------------
  ! Upper limit of the inner y integration.
  !
  ! Change this function to define a different upper boundary y = y(x).
  !---------------------------------------------------------------------
  implicit none

  double precision, intent(in) :: x

  y_upper = sin(x)

end function y_upper


double precision function inner_integral(x)
  !---------------------------------------------------------------------
  ! Evaluate the inner integral
  !
  !             y_upper(x)
  !                 /
  !        F(x) =  |  f(x,y) dy .
  !                 /
  !             y_lower(x)
  !
  ! The resulting value F(x) is supplied to the outer Simpson
  ! integration.
  !
  ! Input:
  !   x - current value of the outer integration variable
  !
  ! Output:
  !   inner_integral - numerical value of the inner y integral
  !---------------------------------------------------------------------
  use nested_data, only: x_current, eps_inner
  implicit none

  double precision :: g, y_lower, y_upper
  double precision, intent(in) :: x
  double precision :: c, d, error_estimate
  integer :: n_used
  external :: g, y_lower, y_upper

  x_current = x

  c = y_lower(x)
  d = y_upper(x)

  call simpson_refine(g, c, d, eps_inner, inner_integral,            &
                      error_estimate, n_used)

end function inner_integral


double precision function g(y)
  !---------------------------------------------------------------------
  ! One-dimensional function seen by the inner Simpson integration.
  !
  ! The current x value is stored in module nested_data.  Thus this
  ! wrapper converts the original f(x,y) into a function of y only.
  !---------------------------------------------------------------------
  use nested_data, only: x_current
  implicit none

  double precision :: integrand
  double precision, intent(in) :: y
  external :: integrand

  g = integrand(x_current,y)

end function g


recursive subroutine simpson_refine(f, a, b, eps, integral,           &
                                    error_estimate, n_used)
  !---------------------------------------------------------------------
  ! One-dimensional Simpson integration with automatic step doubling.
  !
  ! Input:
  !   f   - function to integrate
  !   a   - lower integration limit
  !   b   - upper integration limit
  !   eps - requested tolerance
  !
  ! Output:
  !   integral       - numerical approximation to the integral
  !   error_estimate - |I_(2n)-I_n|/15
  !   n_used         - final number of subintervals
  !
  ! Method:
  !   Composite Simpson estimates are generated with
  !
  !       n = 2, 4, 8, 16, ...
  !
  !   intervals.  For Simpson's rule, the difference between successive
  !   estimates gives the standard error estimate
  !
  !       error approximately = |I_(2n)-I_n|/15.
  !
  !   Refinement stops when this estimate is smaller than eps.
  !
  ! The routine is declared RECURSIVE because, in the present nested
  ! application, the outer Simpson calculation calls inner_integral(x),
  ! which in turn calls this same Simpson routine for the y integration.
  !
  ! Based on code written by Alexander Godunov, October 2009.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: a, b, eps
  double precision, intent(out) :: integral, error_estimate
  integer, intent(out)          :: n_used
  external :: f

  integer, parameter :: nmax = 1048576

  double precision :: sn, s2n, h, x
  integer :: n, i

  if (eps <= 0.0d0) then
    write(*,'(a)') 'simpson_refine: tolerance must be positive.'
    stop
  end if

  if (abs(b-a) <= tiny(1.0d0)) then
    integral = 0.0d0
    error_estimate = 0.0d0
    n_used = 0
    return
  end if

  ! Initial Simpson estimate using two subintervals.
  h = (b-a)/2.0d0

  sn = h*( f(a) + 4.0d0*f(a+h) + f(b) )/3.0d0

  ! Double the number of intervals until the estimated error is small.
  n = 4

  do while (n <= nmax)

    h = (b-a)/dble(n)

    s2n = f(a) + f(b)

    do i = 1, n-1

      x = a + dble(i)*h

      if (mod(i,2) == 0) then
        s2n = s2n + 2.0d0*f(x)
      else
        s2n = s2n + 4.0d0*f(x)
      end if

    end do

    s2n = h*s2n/3.0d0

    error_estimate = abs(s2n-sn)/15.0d0

    if (error_estimate <= eps) then
      integral = s2n
      n_used = n
      return
    end if

    sn = s2n
    n = 2*n

  end do

  ! Maximum refinement reached.
  integral = sn
  n_used = n/2

  write(*,'(a)') 'simpson_refine: maximum number of intervals reached.'

end subroutine simpson_refine