stdlib_quadrature_trapz.fypp Source File

#:include "common.fypp"

submodule (stdlib_quadrature) stdlib_quadrature_trapz
    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
end submodule stdlib_quadrature_trapz