init_piecewise1D Function

    type(piecewise1D) function init_piecewise1D(array) result(res)
        !! initalise the piecewise1D type with an array size (n, 2). Calculates the CDF of this array.
        !> Input array

        use stdlib_quadrature, only: trapz_weights

        real(kind=wp), intent(in) :: array(:, :)
        
        integer :: i, length
        real(kind=wp) :: weights(size(array, 1)), sumer

        if(size(array, 2) /= 2)error stop "Array must be size (n, 2)"

        res%array = array
        length = size(array, 1)
        allocate(res%cdf(length))
        res%cdf = 0.
        ! Generate CDF array from PDF array via Trapezoidal rule
        weights = trapz_weights(array(:, 1))
        sumer = 0.
        do i = 2, length
            sumer = sumer + weights(i)*array(i,2)
            res%cdf(i) = sumer
        end do
         ! normalise
        res%cdf=res%cdf/res%cdf(length)
    end function init_piecewise1D