MODULE quadratures

  IMPLICIT none
  
  integer(kind=8) :: subintervals = 100
  
  INTERFACE
     FUNCTION integrate(ibeg, iend, myfun, p) result(val)
       IMPLICIT none
       ! beginning of integration interval
       real(kind=8), intent(in) :: ibeg
       ! ending of integration interval
       real(kind=8), intent(in) :: iend
       ! function to be integrated
       procedure(funint) :: myfun
       ! polynomial order
       integer(kind=4), intent(in) :: p
       ! result of integration
       real(kind=8) :: val
     END FUNCTION integrate
  END INTERFACE
  
  INTERFACE
     FUNCTION quadrature(qbeg, qend, fun) result(val)
       IMPLICIT none
       real(kind=8), intent(in) :: qbeg, qend
       procedure(funint) :: fun
       real(kind=8) :: val
     END FUNCTION quadrature
  END INTERFACE
  
  INTERFACE
     FUNCTION funint(x) result(y)
       IMPLICIT none
       real(kind=8), intent(in) :: x
       real(kind=8) :: y
     END FUNCTION funint
  END INTERFACE

CONTAINS
  
  FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
    IMPLICIT none
    real(kind=8), intent(in) :: ibeg
    real(kind=8), intent(in) :: iend
    procedure(funint) :: fun
    integer(kind=4), intent(in) :: p
    
    real(kind=8) :: val, subinterval_width, qbeg, qend

    procedure(quadrature), pointer :: quad
    integer(kind=8) :: i

    subinterval_width = (iend - ibeg) / subintervals
    
    SELECT CASE (p)
    CASE (:-1)
       STOP "negative interpolationg polynomial order passed"
    CASE (0)
       quad => rectangle
    CASE (1)
       quad => trapeze
    CASE (2)
       quad => simpson_1_third
    CASE default
       STOP "Newton-Cotes quadratures only implemented for order < 3"
    END SELECT

    ! compute integral using quadrature pointed by quad
    val = 0
    
    DO i = 1, subintervals
       qend = ibeg + i * subinterval_width
       qbeg = ibeg + (i - 1) * subinterval_width
       val = val + quad(qbeg, qend, fun)
    END DO
  END FUNCTION newton_cotes

  FUNCTION rectangle(qbeg, qend, fun) result(val)
    IMPLICIT none
    real(kind=8), intent(in) :: qbeg, qend
    procedure(funint) :: fun
    real(kind=8) :: val

    val = (qend - qbeg) * fun((qend + qbeg) * 0.5)
  END FUNCTION rectangle
  
  FUNCTION trapeze(qbeg, qend, fun) result(val)
    IMPLICIT none
    real(kind=8), intent(in) :: qbeg, qend
    procedure(funint) :: fun
    real(kind=8) :: val

    val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend))
  END FUNCTION trapeze
  
  FUNCTION simpson_1_third(qbeg, qend, fun) result(val)
    IMPLICIT none
    real(kind=8), intent(in) :: qbeg, qend
    procedure(funint) :: fun
    real(kind=8) :: val

    val = (qend - qbeg) * (1/6.0) *                               &
         (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
  END FUNCTION simpson_1_third

END MODULE quadratures