!=======================================================================
! ch07_composite_trapezoidal_simpson.f90
!
! Composite trapezoidal and Simpson rules for numerical integration.
!
! The program evaluates
!
!       integral sin(x) dx,   0 <= x <= pi,
!
! for successively finer uniform grids.  The exact value is 2.
!
! The example illustrates the different convergence rates of the
! composite trapezoidal and Simpson rules as the number of subintervals
! is doubled.
!
! The Simpson routine is based on code written by Alexander Godunov,
! October 2009.  The companion trapezoidal routine and presentation
! were added for the book companion website, 2026.
!=======================================================================

program main
  implicit none

  double precision :: f
  double precision :: a, b, exact
  double precision :: trap, simp
  double precision, parameter :: pi = 3.14159265358979323846d0
  integer :: n, i
  external :: f

  a = 0.0d0
  b = pi
  exact = 2.0d0

  write(*,'(a)')        'Composite trapezoidal and Simpson rules'
  write(*,'(a)')        'Integral of sin(x) from 0 to pi'
  write(*,'(a,es15.7)') 'Exact value          = ', exact
  write(*,*)
  write(*,'(a)') '        n       Trapezoidal          Simpson'
  write(*,'(a)') '------------------------------------------------'

  n = 2

  do i = 1, 16

    call composite_trapezoidal(f, a, b, n, trap)
    call composite_simpson(f, a, b, n, simp)

    write(*,'(i9,2es19.7)') n, trap, simp

    n = 2*n

  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 composite_trapezoidal(f, a, b, n, integral)
  !---------------------------------------------------------------------
  ! Composite trapezoidal rule for integration of f(x) on [a,b].
  !
  ! Input:
  !   f - function to integrate
  !   a - lower integration limit
  !   b - upper integration limit
  !   n - number of equal subintervals
  !
  ! Output:
  !   integral - numerical approximation to the integral
  !
  ! Method:
  !   The interval [a,b] is divided into n equal subintervals of width
  !
  !       h = (b-a)/n .
  !
  !   The composite trapezoidal rule is
  !
  !       integral = h [ f(a)/2 + sum f(x_i) + f(b)/2 ],
  !
  !   where x_i = a + i h and i = 1,...,n-1.
  !---------------------------------------------------------------------
  implicit none

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

  double precision :: h, x, sum
  integer :: i
  external :: f

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

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

  sum = 0.5d0*(f(a) + f(b))

  do i = 1, n-1
    x = a + dble(i)*h
    sum = sum + f(x)
  end do

  integral = h*sum

end subroutine composite_trapezoidal


subroutine composite_simpson(f, a, b, n, integral)
  !---------------------------------------------------------------------
  ! Composite Simpson rule for integration of f(x) on [a,b].
  !
  ! Input:
  !   f - function to integrate
  !   a - lower integration limit
  !   b - upper integration limit
  !   n - number of equal subintervals; n must be even
  !
  ! Output:
  !   integral - numerical approximation to the integral
  !
  ! Method:
  !   The interval [a,b] is divided into n equal subintervals of width
  !
  !       h = (b-a)/n .
  !
  !   The composite Simpson rule is
  !
  !       integral = h/3 [ f(a) + f(b)
  !                        + 4 sum f(x_i), i odd
  !                        + 2 sum f(x_i), i even ] .
  !
  !   The number of subintervals must be even.
  !
  ! Based on code 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

  double precision :: h, x
  double precision :: sum_odd, sum_even
  integer :: i
  external :: f

  if (n < 2 .or. mod(n,2) /= 0) then
    write(*,'(a)') 'composite_simpson: n must be a positive even integer.'
    stop
  end if

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

  sum_odd = 0.0d0
  do i = 1, n-1, 2
    x = a + dble(i)*h
    sum_odd = sum_odd + f(x)
  end do

  sum_even = 0.0d0
  do i = 2, n-2, 2
    x = a + dble(i)*h
    sum_even = sum_even + f(x)
  end do

  integral = h*(f(a) + f(b) + 4.0d0*sum_odd + 2.0d0*sum_even)/3.0d0

end subroutine composite_simpson