!=======================================================================
! ch07_principal_value.f90
!
! Cauchy principal value integral using analytical subtraction.
!
! The program evaluates
!
!                 1
!                 /
!                |       cos(x)
!      PV        |     --------- dx
!                |      x - 0.2
!                /
!                -1
!
! The singularity at x0 = 0.2 is removed analytically by writing
!
!        f(x)        f(x)-f(x0)        f(x0)
!      ---------  =  -----------  +  --------- .
!       x-x0           x-x0            x-x0
!
! Therefore
!
!                                      b-x0
! PV integral = integral [f(x)-f(x0)]/(x-x0) dx
!               + f(x0) log(-------------) .
!                                      x0-a
!
! The remaining numerical integral is regular.  It is evaluated on
! [a,x0] and [x0,b] with Gauss-Legendre quadrature.  Gauss-Legendre
! nodes do not include interval endpoints, so x = x0 is never sampled.
!
! This example follows Section 7.5, "Principal value integrals," of
! Computational Physics: From Equations to Simulation.
!
! Alexander Godunov
! Companion website, 2026
!=======================================================================

program main
  implicit none

  double precision :: f
  double precision :: a, b, x0, result, reference
  integer :: n, i
  integer, dimension(4) :: nvalues
  external :: f

  a  = -1.0d0
  b  =  1.0d0
  x0 =  0.2d0

  reference = -0.5912784964342436d0

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

  write(*,'(a)')        'Cauchy principal value integral'
  write(*,'(a)')        'PV integral of cos(x)/(x-0.2) from -1 to 1'
  write(*,'(a,es15.7)') 'Reference value      = ', reference
  write(*,*)
  write(*,'(a)') '        n       Principal value'
  write(*,'(a)') '--------------------------------'

  do i = 1, size(nvalues)

    n = nvalues(i)

    call principal_value(f, a, b, x0, n, result)

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

  end do

end program main


double precision function f(x)
  !---------------------------------------------------------------------
  ! Numerator function in the principal-value integral
  !
  !                 b
  !                 /
  !                |       f(x)
  !      PV        |     -------- dx .
  !                |      x - x0
  !                /
  !                 a
  !
  ! To use this program for another problem of the same form, replace
  ! only the line defining f(x) below and set a, b, and x0 in the main
  ! program.  The principal-value and quadrature routines do not need
  ! to be changed.
  !
  ! Input:
  !   x - integration variable
  !
  ! Output:
  !   f - value of the numerator function
  !---------------------------------------------------------------------
  implicit none

  double precision, intent(in) :: x

  f = cos(x)

end function f


subroutine principal_value(f, a, b, x0, n, result)
  !---------------------------------------------------------------------
  ! Cauchy principal value of
  !
  !                 b
  !                 /
  !                |       f(x)
  !      PV        |     -------- dx
  !                |      x - x0
  !                /
  !                 a
  !
  ! using analytical subtraction of the singularity.
  !
  ! Input:
  !   f  - numerator function
  !   a  - lower integration limit
  !   b  - upper integration limit
  !   x0 - location of the pole, a < x0 < b
  !   n  - number of Gauss-Legendre points on each subinterval
  !
  ! Output:
  !   result - numerical value of the principal-value integral
  !
  ! Method:
  !
  !   Write
  !
  !        f(x)        f(x)-f(x0)        f(x0)
  !      ---------  =  -----------  +  --------- .
  !       x-x0           x-x0            x-x0
  !
  !   The first term is regular at x0, while the principal value of
  !   the second term is known analytically:
  !
  !                                   b-x0
  !      PV integral dx/(x-x0) = log(-------) .
  !                                   x0-a
  !
  !   The regular part is integrated numerically on [a,x0] and
  !   [x0,b].  Splitting the interval keeps x0 at an endpoint, and
  !   Gauss-Legendre quadrature does not evaluate interval endpoints.
  !---------------------------------------------------------------------
  implicit none

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

  double precision :: f0, left, right, singular_part

  if (x0 <= a .or. x0 >= b) then
    write(*,'(a)') 'principal_value: x0 must lie inside (a,b).'
    stop
  end if

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

  f0 = f(x0)

  call gauss_legendre_regular(f, a,  x0, x0, f0, n, left)
  call gauss_legendre_regular(f, x0, b,  x0, f0, n, right)

  singular_part = f0*log((b-x0)/(x0-a))

  result = left + right + singular_part

end subroutine principal_value


subroutine gauss_legendre_regular(f, a, b, x0, f0, n, integral)
  !---------------------------------------------------------------------
  ! Gauss-Legendre quadrature for the regularized function
  !
  !             f(x)-f(x0)
  !             -----------
  !                x-x0
  !
  ! on a finite interval [a,b].
  !
  ! Input:
  !   f  - numerator function
  !   a  - lower integration limit
  !   b  - upper integration limit
  !   x0 - location of the original pole
  !   f0 - f(x0)
  !   n  - number of Gauss-Legendre points
  !
  ! Output:
  !   integral - numerical integral of the regularized function
  !
  ! Method:
  !   The Gauss-Legendre nodes are roots of the Legendre polynomial
  !   P_n.  They are found by Newton iteration.  The corresponding
  !   weights are calculated from P'_n.
  !
  !   The standard nodes on [-1,1] are mapped to [a,b].  In the present
  !   application x0 is an endpoint of each subinterval, and therefore
  !   is never sampled by Gauss-Legendre quadrature.
  !
  ! Based on the original Gauss quadrature routines written by
  ! Alexander Godunov, October 2009.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: a, b, x0, f0
  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
  double precision :: gleft, gright
  integer :: i, j, m, iter

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

  m = (n+1)/2

  do i = 1, m

    ! Initial approximation to a root of P_n.
    z = cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

    ! Newton iteration for the root.
    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_regular: Newton iteration failed.'
      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

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

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

    gleft  = (f(xleft) -f0)/(xleft -x0)
    gright = (f(xright)-f0)/(xright-x0)

    if (abs(z) <= tol) then
      integral = integral + weight*(f(midpoint)-f0)/(midpoint-x0)
    else
      integral = integral + weight*(gleft+gright)
    end if

  end do

  integral = halfwidth*integral

end subroutine gauss_legendre_regular