! 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 H05M2
    use H05M1  ! Defines function F and parameter HI
    implicit none
    private
    public :: Quadrature, Quad, Simpson
    type, public :: Result_Type
      real(kind = HI) :: X, Err_Loc
      logical :: Error
    end type Result_Type
  contains

    function Quadrature( A, B, Tol, M_Int ) result ( Quadrature_R )
      real(kind = HI), intent(in) :: A, B, Tol, M_Int
      type(Result_Type) :: Quadrature_R
      real(kind = HI), parameter :: Magic_No = 10.0   ! See Burden & Faires, p. 195.
  ! start function Quadrature
      Quadrature_R = Quad( A, B, Simpson( A, B ), Tol * Magic_No, M_Int )
      return
    end function Quadrature

    recursive function Quad( A, B, I, Tol, M_Int ) result( Quad_R )
      real(kind = HI), intent (in) :: A, B, I, Tol, M_Int
      type(Result_Type) :: Quad_R
      type(Result_Type) :: R1, R2
      real(kind = HI) :: Y, Z, M
  ! start function Quad
      M = (A + B) / 2.0_HI
      Quad_R = Result_Type( I, M, .false. )
      Y = Simpson( A, M )
      Z = Simpson( M, B )
      if (abs( Y + Z - Quad_R % X ) > Tol) then
        if (abs( B - A ) < M_Int) then ! Not OK but interval is too small
          Quad_R % Error = .true.
        else ! Not OK but keep trying
          R1 = Quad( A, M, Y, Tol / 2.0_HI, M_Int )
          R2 = Quad( M, B, Z, Tol / 2.0_HI, M_Int )
          Quad_R % X = R1 % X + R2 % X
          Quad_R % Error = R1 % Error .or. R2 % Error
          if (Quad_R % Error) then
            Quad_R % Err_Loc = merge(R1 % Err_Loc, R2 % Err_Loc, R1 % Error)
          end if
        end if  
      else ! OK
        Quad_R % X = Y + Z
      end if
      return
    end function Quad

    function Simpson( A, B ) result( Simpson_R )
      real (kind = HI), intent(in) :: A, B
      real (kind = HI) :: Simpson_R
      real (kind = HI) :: M
  ! start function Simpson
      M = (A + B) / 2.0_HI
      Simpson_R = (B - A) * (F( A ) + 4.0_HI * F( M ) + F( B )) / 6.0_HI
      return
    end function Simpson

  end module H05M2

