!=======================================================================
! ch07_gauss_legendre.f90
!
! Gauss-Legendre quadrature for numerical integration.
!
! The program evaluates
!
!       integral sin(x) dx,   0 <= x <= pi,
!
! using several n-point Gauss-Legendre rules.  The exact value is 2.
!
! Unlike a version with tabulated nodes and weights for only selected
! values of n, this program computes the Gauss-Legendre nodes and
! weights numerically and can therefore be used for any positive n.
!
! The roots of the Legendre polynomial P_n(x) are found by Newton
! iteration.  Symmetry is used so that only half of the roots need to
! be computed explicitly.
!
! Based on the original 8-point and 16-point Gauss quadrature programs
! written by Alexander Godunov, October 2009.
! Revised and generalized for the companion website, 2026.
!=======================================================================

program main
  implicit none

  double precision :: f
  double precision :: a, b, exact, integral
  double precision, parameter :: pi = 3.14159265358979323846d0
  integer :: n, i
  integer, dimension(6) :: nvalues
  external :: f

  a = 0.0d0
  b = pi
  exact = 2.0d0

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

  write(*,'(a)')        'Gauss-Legendre quadrature'
  write(*,'(a)')        'Integral of sin(x) from 0 to pi'
  write(*,'(a,es15.7)') 'Exact value          = ', exact
  write(*,*)
  write(*,'(a)') '        n       Gauss-Legendre'
  write(*,'(a)') '--------------------------------'

  do i = 1, size(nvalues)

    n = nvalues(i)
    call gauss_legendre(f, a, b, n, integral)

    write(*,'(i9,es21.7)') n, integral

  end do

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 sin(x)
  !---------------------------------------------------------------------
  implicit none

  double precision, intent(in) :: x

  f = sin(x)

end function f


subroutine gauss_legendre(f, a, b, n, integral)
  !---------------------------------------------------------------------
  ! n-point Gauss-Legendre quadrature for integration of f(x) on [a,b].
  !
  ! Input:
  !   f - function to integrate
  !   a - lower integration limit
  !   b - upper integration limit
  !   n - number of Gauss-Legendre points
  !
  ! Output:
  !   integral - numerical approximation to the integral
  !
  ! Method:
  !   Gauss-Legendre quadrature is defined on the standard interval
  !   [-1,1].  Its nodes are the roots of the Legendre polynomial P_n.
  !   The roots are found here by Newton iteration.
  !
  !   The Legendre polynomial is evaluated by 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 .
  !
  !   At a root z of P_n, the derivative is
  !
  !       P'_n(z) = n [ z P_n(z) - P_{n-1}(z) ] / (z^2 - 1),
  !
  !   and the corresponding quadrature weight is
  !
  !       w = 2 / [ (1-z^2) (P'_n(z))^2 ] .
  !
  !   The nodes and weights are symmetric about zero, so only half of
  !   them are computed explicitly.  The transformation
  !
  !       x = (b-a)z/2 + (a+b)/2
  !
  !   maps each node z from [-1,1] to the interval [a,b].
  !
  !   For sufficiently smooth functions, increasing n produces very
  !   rapid convergence compared with equally spaced Newton-Cotes rules.
  !
  ! Based on the original fixed 8-point and 16-point routines written
  ! by Alexander Godunov, October 2009.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: a, b
  integer, intent(in)           :: n
  double precision, intent(out) :: integral
  external :: f

  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 :: weight
  double precision :: midpoint, halfwidth
  double precision :: xleft, xright
  integer :: i, j, m, iter

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

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

  ! Because the roots are symmetric, only the positive half is needed.
  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(z).
    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: Newton iteration did not converge.'
      stop
    end if

    ! Re-evaluate P_n and P_{n-1} at the converged root so that the
    ! derivative and weight correspond to the final value of z.
    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

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

    xleft  = midpoint - halfwidth*z
    xright = midpoint + halfwidth*z

    ! For odd n the central root is z = 0 and must be counted once.
    if (abs(z) <= tol) then
      integral = integral + weight*f(midpoint)
    else
      integral = integral + weight*(f(xleft) + f(xright))
    end if

  end do

  integral = halfwidth*integral

end subroutine gauss_legendre