
MODULE integration_mod

  IMPLICIT NONE

!-----------------------------------------------------------------------------
! Description:
!   This module numerically integrates mathematical functions
! This module is added here in order to calculate integrals
! http://rosettacode.org/wiki/Numerical_integration#Fortran
!
! Current Code Owner: Renato Braghiere
!
! Code Description:
!   Language: Fortran 90.
!   This code is written to JULES coding standards v1.
!-----------------------------------------------------------------------------

!-----------------------------------------------------------------------------
! Module constants
!-----------------------------------------------------------------------------
 
CONTAINS
 
  ! function, lower limit, upper limit, steps, method
  FUNCTION integrate(f, a, b, in, method)
    real :: integrate
    real, intent(in) :: a, b
    integer, optional, intent(in) :: in
    character(len=*), intent(in), optional :: method
    interface
     elemental function f(ra)
    
         real :: f
         real, intent(in) :: ra
       end function f
    end interface
 
    integer :: n, i, m
    real :: h
    real, dimension(:), allocatable :: xpoints
    real, dimension(:), target, allocatable :: fpoints
    real, dimension(:), pointer :: fleft, fmid, fright
 
    if ( present(in) ) then
       n = in
    else
       n = 20
    end if
 
    if ( present(method) ) then
       select case (method)
       case ('leftrect')
          m = 1
       case ('midrect')
          m = 2
       case ('rightrect')
          m = 3
       case ( 'trapezoid' )
          m = 4
       case default
          m = 0
       end select
    else
       m = 0
    end if
 
    h = (b - a) / n
 
    allocate(xpoints(0:2*n), fpoints(0:2*n))
 
    xpoints = (/ (a + h*i/2, i = 0,2*n) /)
 
    fpoints = f(xpoints)
    fleft  => fpoints(0 : 2*n-2 : 2)
    fmid   => fpoints(1 : 2*n-1 : 2)
    fright => fpoints(2 : 2*n   : 2)
 
    select case (m)
    case (0) ! simpson
       integrate = h / 6.0 * sum(fleft + fright + 4.0*fmid)
    case (1) ! leftrect
       integrate = h * sum(fleft)
    case (2) ! midrect
       integrate = h * sum(fmid)
    case (3) ! rightrect
       integrate = h * sum(fright)
    case (4) ! trapezoid
       integrate = h * sum(fleft + fright) / 2
    end select
 
    deallocate(xpoints, fpoints)
  END FUNCTION integrate
 
END MODULE integration_mod



