! Program Example from "Essential Fortran 90 & 95" by Loren P. Meissner
! Copyright 1996. Copying for sale requires permission from the author.
! Otherwise, distribution is permitted if these three lines are included.

  module I07M1 !  Perform quadrature by Simpsons rule.
    implicit none
    public :: Quadrature
  contains
    function Quadrature( Integrand, A, B, N_In ) result( Simpson )
      interface
        function Integrand( X ) result( Y )
          real, intent (in) :: X
          real :: Y
        end function Integrand
      end interface
      real, intent(in) :: A, B
      integer, intent (in) :: N_In
      real :: Simpson
      real :: Evens, Odds, X, Delta_X
      integer :: I, N
  ! start function Quadrature
      N = N_In + modulo( N_In, 2 )   ! Add 1 to N if it is not even.
      Delta_X = (B - A) / real( N )
      Evens = 0.0
      Odds = 0.0
      do I = 1, N - 1
        X = A + real( I ) * Delta_X
        if (modulo( I, 2 ) == 0) then
          Evens = Evens + Integrand( X )
        else
          Odds = Odds + Integrand( X )
        end if
      end do
      Simpson = Delta_X * &
        (Integrand( A ) + Integrand( B ) + 4.0 * Odds + 2.0 * Evens) / 3.0
      return
    end function Quadrature
  end module I07M1

