integrate1D Module Function

    module function integrate1D(fun, xbnds) result(result)
        ! This routine uses gauss-legendre quadrature to integrate a 1D
        ! line function

        ! Input/Output dummy variables
        real(wp), intent(in)    :: xbnds(2)
        real(wp)                :: result
        procedure(fun1d_interf) :: fun

        ! Local variables
        integer, parameter    :: N_start = 2
        integer               :: N
        real(wp)              :: result_old, error
        real(wp), parameter   :: eps = epsilon(0e0)
        real(wp), allocatable :: x(:), w(:), y(:)

        N = N_start
        result = 0.0_wp
        error = 1.0_wp

        do
            result_old = result

            allocate(x(N), w(N), y(N))

            call gaussquad(N, x, w)

            y = fun(x)
            result = dot_product(y, w)

            deallocate(x, w, y)

            ! error = norm2([result, result_old])
            error = sqrt((result-result_old)**2.0_wp)
            ! error = abs((result-result_old)/(result_old+eps))
            ! print*, N, result, error
            ! print*, N, result, result_old, error, eps

            ! Check if error is acceptable, as well as whether the loop
            ! was gone through at least twice (N = 3, 4, 5...)
            if ( N > N_start .and. &
                norm2( [ error ] ) <= eps ) then
                ! print'(a,i3,a,e13.5)',  'Fun integrated in ', N, &
                !                         ' iterations. Error = ', error
                exit
            else
                N = N + 1
            endif
        enddo
        return
    end function integrate1D