!_______________________________________________________________________
! This module contains numerical integration routines
! from the module adapt_quad by Allan Miller
! Aleksandar Donev, November 2000
!_______________________________________________________________________
MODULE Numerical_Integration
      USE Adapt_Quad
      PRIVATE
!
      PUBLIC :: AdaptiveIntegral
      INTEGER, PARAMETER, PUBLIC :: wp=KIND(0.0D0)
      INTEGER, SAVE, PUBLIC :: error_flag = 0
!
CONTAINS
  !___________________________________
  !            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A
  !            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
  !            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  !            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  !___________________________________
      FUNCTION AdaptiveIntegral (f, a, b, error_estimate, relative_error, absolute_error) RESULT (integral)
!
         ! We must supply an interface for the function
         INTERFACE
            FUNCTION f (x) RESULT (f_x)! This would be a pure function in Fortran 90
               REAL (KIND=KIND(0.0D0)), INTENT (IN) :: x
               REAL (KIND=KIND(0.0D0)) :: f_x
            END FUNCTION f
         END INTERFACE
         REAL (KIND=wp), INTENT (IN) :: a, b
         REAL (KIND=wp), INTENT (OUT) :: error_estimate
         REAL (KIND=wp), INTENT (IN), OPTIONAL :: relative_error, absolute_error
         REAL (KIND=wp) :: integral
!
         REAL (KIND=wp) :: abserr, relerr
         INTEGER, PARAMETER :: n_intervals = 100
         INTEGER :: error_status = 0, last = 0
!
         IF (present(absolute_error)) THEN
            abserr = absolute_error
         ELSE
            abserr = 1.0E-3
         END IF
!
         IF (present(relative_error)) THEN
            relerr = relative_error
         ELSE
            relerr = 1.0E-3
         END IF
!
         CALL qxgs (f, a, b, abserr, relerr, integral, error_estimate, error_status, n_intervals, last)
         ! I only raise the error flag, it is the users duty to lower it if needed
         IF (error_status /= 0) error_flag = error_status
!
      END FUNCTION AdaptiveIntegral
!
END MODULE Numerical_Integration
!_______________________________________________________________________