## Source Code

#:include "common.fypp"

use stdlib_error, only: check
implicit none

contains

#:for KIND in REAL_KINDS

pure module function trapz_dx_${KIND}$(y, dx) result(integral)
real(${KIND}$), dimension(:), intent(in) :: y
real(${KIND}$), intent(in) :: dx
real(${KIND}$) :: integral

integer :: n

n = size(y)

select case (n)
case (0:1)
integral = 0.0_${KIND}$
case (2)
integral = 0.5_${KIND}$*dx*(y(1) + y(2))
case default
integral = dx*(sum(y(2:n-1)) + 0.5_${KIND}$*(y(1) + y(n)))
end select
end function trapz_dx_${KIND}$

#:endfor
#:for KIND in REAL_KINDS

module function trapz_x_${KIND}$(y, x) result(integral)
real(${KIND}$), dimension(:), intent(in) :: y
real(${KIND}$), dimension(:), intent(in) :: x
real(${KIND}$) :: integral

integer :: i
integer :: n

n = size(y)
call check(size(x) == n, "trapz: Arguments x and y must be the same size.")

select case (n)
case (0:1)
integral = 0.0_${KIND}$
case (2)
integral = 0.5_${KIND}$*(x(2) - x(1))*(y(1) + y(2))
case default
integral = 0.0_${KIND}$
do i = 2, n
integral = integral + (x(i) - x(i-1))*(y(i) + y(i-1))
end do
integral = 0.5_${KIND}$*integral
end select
end function trapz_x_${KIND}$

#:endfor
#:for KIND in REAL_KINDS

pure module function trapz_weights_${KIND}$(x) result(w)
real(${KIND}$), dimension(:), intent(in) :: x
real(${KIND}$), dimension(size(x)) :: w

integer :: i
integer :: n

n = size(x)

select case (n)
case (0)
! no action needed
case (1)
w(1) = 0.0_${KIND}$
case (2)
w = 0.5_${KIND}$*(x(2) - x(1))
case default
w(1) = 0.5_${KIND}$*(x(2) - x(1))
w(n) = 0.5_${KIND}$*(x(n) - x(n-1))
do i = 2, size(x)-1
w(i) = 0.5_${KIND}$*(x(i+1) - x(i-1))
end do
end select
end function trapz_weights_${KIND}$

#:endfor