integration.f90 Source File

!> Utility module containing generic integration routines
module integration
  implicit none
  private

  public :: trapezoidal_integration
  public :: rectangular_integration

  interface trapezoidal_integration
     module procedure real_trapezoidal_integration
     module procedure complex_trapezoidal_integration
  end interface trapezoidal_integration

  interface rectangular_integration
     module procedure real_rectangular_integration
     module procedure complex_rectangular_integration
  end interface rectangular_integration

  interface local_sum
     module procedure real_local_sum
     module procedure complex_local_sum
  end interface local_sum

  logical, parameter :: use_compensated_sum = .true.
contains
  !> A wrapper routine to perform a sum. Used to provide a single
  !> point at which we choose a summation method
  real pure function real_local_sum(array) result(total)
    use summation, only: klein_sum
    implicit none
    real, dimension(:), intent(in) :: array

    if (use_compensated_sum) then
       total = klein_sum(array)
    else
       total = sum(array)
    end if
  end function real_local_sum

  !> A wrapper routine to perform a sum. Used to provide a single
  !> point at which we choose a summation method
  complex pure function complex_local_sum(array) result(total)
    use summation, only: klein_sum
    implicit none
    complex, dimension(:), intent(in) :: array

    if (use_compensated_sum) then
       total = cmplx(klein_sum(real(array)),klein_sum(aimag(array)))
    else
       total = sum(array)
    end if
  end function complex_local_sum

  !> Apply the trapezium rule to integrate y over the domain of x.
  real pure function real_trapezoidal_integration(x, y) result(integral)
    implicit none
    real, dimension(:), intent(in) :: x, y
    integer :: n
    n = size(x)
    integral = local_sum((x(2:)-x(:n-1))*(y(:n-1)+y(2:))) * 0.5
  end function real_trapezoidal_integration

  !> Apply the trapezium rule to integrate y over the domain of x.
  complex pure function complex_trapezoidal_integration(x, y) result(integral)
    implicit none
    real, dimension(:), intent(in) :: x
    complex, dimension(:), intent(in) :: y
    integer :: n
    n = size(x)
    integral = local_sum((x(2:)-x(:n-1))*(y(:n-1)+y(2:))) * 0.5
  end function complex_trapezoidal_integration

  !> Apply the (left) rectangular rule to integrate y over the domain of x.
  real pure function real_rectangular_integration(x, y) result(integral)
    implicit none
    real, dimension(:), intent(in) :: x, y
    integer :: n
    n = size(x)
    integral = local_sum((x(2:)-x(:n-1))*y(:n-1))
  end function real_rectangular_integration

  !> Apply the (left) rectangular rule to integrate y over the domain of x.
  complex pure function complex_rectangular_integration(x, y) result(integral)
    implicit none
    real, dimension(:), intent(in) :: x
    complex, dimension(:), intent(in) :: y
    integer :: n
    n = size(x)
    integral = local_sum((x(2:)-x(:n-1))*y(:n-1))
  end function complex_rectangular_integration

end module integration