!=======================================================================
! ch07_gauss_legendre_2d.f90
!
! Two-dimensional integration using a Gauss-Legendre product rule.
!
! The program evaluates
!
!          1   1
!          /   /
!         |   |        1
!     I = |   |   -------------  dy dx
!         |   |   1 + x^2 + y^2
!         /   /
!          0   0
!
! This example follows the discussion of product rules in Section 7.7,
! "Multidimensional integration," of
! Computational Physics: From Equations to Simulation.
!
! An n-point Gauss-Legendre rule is used in each coordinate direction.
! The resulting two-dimensional product rule uses n^2 function
! evaluations:
!
!          I approximately = sum_i sum_j wi wj f(xi,yj).
!
! The example therefore illustrates both the rapid convergence of
! Gaussian quadrature for a smooth low-dimensional problem and the
! growth in computational cost as the number of dimensions increases.
!
! Alexander Godunov
! Companion website, 2026
!=======================================================================

program main
  implicit none

  double precision :: f
  double precision :: ax, bx, ay, by
  double precision :: integral, reference, error
  integer :: n, nfun, i
  integer, dimension(5) :: nvalues
  external :: f

  ax = 0.0d0
  bx = 1.0d0
  ay = 0.0d0
  by = 1.0d0

  ! High-accuracy reference value for the test integral.
  ! If f(x,y) or the integration limits are changed, this value should
  ! also be replaced or omitted.
  reference = 6.3951035187031102d-01

  nvalues = (/ 2, 4, 8, 16, 32 /)

  write(*,'(a)')        'Two-dimensional Gauss-Legendre integration'
  write(*,'(a)')        'Integral of 1/(1+x^2+y^2) over [0,1] x [0,1]'
  write(*,'(a,es15.7)') 'Reference value      = ', reference
  write(*,*)
  write(*,'(a)') '        n     Function calls          Integral       Absolute error'
  write(*,'(a)') '-------------------------------------------------------------------'

  do i = 1, size(nvalues)

    n = nvalues(i)

    call gauss_legendre_2d(f, ax, bx, ay, by, n, integral, nfun)

    error = abs(integral-reference)

    write(*,'(i9,i19,2es20.7)') n, nfun, integral, error

  end do

end program main


double precision function f(x,y)
  !---------------------------------------------------------------------
  ! Function to be integrated over the rectangular domain
  !
  !          bx  by
  !          /   /
  !         |   |
  !     I = |   |  f(x,y) dy dx .
  !         |   |
  !         /   /
  !          ax  ay
  !
  ! To use this program for another two-dimensional problem over a
  ! rectangular domain, replace only the line defining f(x,y) below
  ! and change the limits in the main program.
  !
  ! Input:
  !   x, y - integration variables
  !
  ! Output:
  !   f    - value of the two-dimensional integrand
  !---------------------------------------------------------------------
  implicit none

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

  f = 1.0d0/(1.0d0 + x*x + y*y)

end function f


subroutine gauss_legendre_2d(f, ax, bx, ay, by, n, integral, nfun)
  !---------------------------------------------------------------------
  ! Two-dimensional Gauss-Legendre product rule on a rectangle.
  !
  ! Input:
  !   f      - function f(x,y) to integrate
  !   ax,bx  - lower and upper limits in x
  !   ay,by  - lower and upper limits in y
  !   n      - number of Gauss-Legendre points in each direction
  !
  ! Output:
  !   integral - numerical approximation to the double integral
  !   nfun     - number of function evaluations
  !
  ! Method:
  !   An n-point Gauss-Legendre rule is generated independently in
  !   the x and y directions.  The two-dimensional quadrature is the
  !   tensor product of the one-dimensional rules:
  !
  !       integral = sum_i sum_j wx(i) wy(j) f(x(i),y(j)).
  !
  !   Thus n points in each direction require n^2 evaluations.
  !   In d dimensions, the corresponding product rule requires n^d
  !   evaluations, illustrating the curse of dimensionality.
  !
  ! This product-rule approach is most useful for smooth functions in
  ! low dimensions, particularly d <= 3.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: ax, bx, ay, by
  integer, intent(in)           :: n
  double precision, intent(out) :: integral
  integer, intent(out)          :: nfun
  external :: f

  double precision :: xnode(n), xweight(n)
  double precision :: ynode(n), yweight(n)
  integer :: i, j

  if (n < 1) then
    write(*,'(a)') 'gauss_legendre_2d: n must be positive.'
    stop
  end if

  call gauss_legendre_nodes_weights(n, ax, bx, xnode, xweight)
  call gauss_legendre_nodes_weights(n, ay, by, ynode, yweight)

  integral = 0.0d0
  nfun = 0

  do i = 1, n
    do j = 1, n

      integral = integral                                   &
               + xweight(i)*yweight(j)*f(xnode(i),ynode(j))

      nfun = nfun + 1

    end do
  end do

end subroutine gauss_legendre_2d


subroutine gauss_legendre_nodes_weights(n, a, b, node, weight)
  !---------------------------------------------------------------------
  ! Generate n Gauss-Legendre nodes and weights on the interval [a,b].
  !
  ! Input:
  !   n      - number of quadrature points
  !   a, b   - lower and upper integration limits
  !
  ! Output:
  !   node   - Gauss-Legendre nodes mapped to [a,b]
  !   weight - corresponding weights on [a,b]
  !
  ! Method:
  !   The Gauss-Legendre nodes are the roots of the Legendre polynomial
  !   P_n(x).  The roots are found by Newton iteration.
  !
  !   The Legendre polynomial is evaluated using the recurrence
  !
  !       P_0(x) = 1,
  !       P_1(x) = x,
  !
  !       P_j(x) = [(2j-1)x P_(j-1)(x) - (j-1)P_(j-2)(x)]/j.
  !
  !   The nodes and weights are first obtained on [-1,1] and then
  !   mapped to the requested interval [a,b].
  !
  ! Based on the one-dimensional Gauss quadrature routines written by
  ! Alexander Godunov, October 2009.
  !---------------------------------------------------------------------
  implicit none

  integer, intent(in) :: n
  double precision, intent(in) :: a, b
  double precision, intent(out) :: node(n), weight(n)

  double precision, parameter :: pi = 3.14159265358979323846d0
  double precision, parameter :: tol = 1.0d-14
  integer, parameter :: max_iter = 100

  double precision :: z, zold
  double precision :: p0, p1, p2, dp
  double precision :: w
  double precision :: midpoint, halfwidth
  integer :: i, j, m, iter

  if (n < 1) then
    write(*,'(a)') 'gauss_legendre_nodes_weights: n must be positive.'
    stop
  end if

  midpoint  = (a+b)/2.0d0
  halfwidth = (b-a)/2.0d0

  ! Because the roots are symmetric, only half are computed directly.
  m = (n+1)/2

  do i = 1, m

    ! Initial approximation to the i-th positive root.
    z = cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

    ! Newton iteration for a root of P_n.
    do iter = 1, max_iter

      p0 = 1.0d0
      p1 = z

      if (n == 1) then
        p2 = p1
      else
        do j = 2, n
          p2 = ((2.0d0*dble(j)-1.0d0)*z*p1              &
               - (dble(j)-1.0d0)*p0)/dble(j)
          p0 = p1
          p1 = p2
        end do
      end if

      if (n == 1) then
        dp = 1.0d0
      else
        dp = dble(n)*(z*p1-p0)/(z*z-1.0d0)
      end if

      zold = z
      z = zold - p1/dp

      if (abs(z-zold) <= tol) exit

    end do

    if (iter > max_iter) then
      write(*,'(a)') 'Gauss-Legendre root iteration did not converge.'
      stop
    end if

    ! Re-evaluate P_n and P_(n-1) at the converged root.
    p0 = 1.0d0
    p1 = z

    if (n == 1) then
      dp = 1.0d0
    else
      do j = 2, n
        p2 = ((2.0d0*dble(j)-1.0d0)*z*p1              &
             - (dble(j)-1.0d0)*p0)/dble(j)
        p0 = p1
        p1 = p2
      end do

      dp = dble(n)*(z*p1-p0)/(z*z-1.0d0)
    end if

    w = 2.0d0/((1.0d0-z*z)*dp*dp)

    ! Symmetric nodes and weights mapped from [-1,1] to [a,b].
    node(i)       = midpoint - halfwidth*z
    node(n+1-i)   = midpoint + halfwidth*z

    weight(i)     = halfwidth*w
    weight(n+1-i) = halfwidth*w

  end do

end subroutine gauss_legendre_nodes_weights