integrate_basis_1d_Ie Function

    function integrate_basis_1d_Ie(N, ii, jj, dii, djj, xy) result(integral)
        integer,  intent(in)                    :: N, ii, jj, dii, djj
        real(wp), intent(in), dimension(:)      :: xy
        real(wp)                                :: integral

        integer                                 :: kk
        integer                                 :: order, basis_num1, basis_num2
        real(wp), dimension(N+1,N+1)            :: Vinv

        order = N
        basis_num1 = ii
        basis_num2 = jj

        Vinv = getAlpha(order+1)

        ! Check to make sure xy is an array of size (N+1)
        if ( size(xy) /=  order+1 ) then
            write(*,*) 'The shape of `xy` is outside the acceptable range for a'
            write(*,*) '1D basis function.'
            write(*,*) 'Shape(xy) should be [', 2, order+1, '], not ', shape(xy)
        endif

        ! call integrate(local_wrapper, [-1.0_wp, 1.0_wp], integral)
        integral = integrate(local_wrapper, [-1._wp, 1._wp])

        return
    contains
        pure function local_wrapper(s) result(y)
            real(wp), intent(in)  :: s(:)
            real(wp), allocatable :: y(:)

            real(wp), allocatable :: J(:)

            allocate(J, mold=s)
            allocate(y, mold=s)
            call XorJ(s, 1, J)

            y = 1.0_wp

            ! Here we have to be careful because J is not always needed in the
            ! first two function calls. Instead of using if statements, we can
            ! use an exponent so that when dx_ == 0, J is 1
            y = y * basis_1D(s, Vinv(:, basis_num1), dii) / (J**dble(dii))
            y = y * basis_1D(s, Vinv(:, basis_num2), djj) / (J**dble(djj))
            y = y * J

            deallocate(J)

            return
        end function local_wrapper

        pure subroutine XorJ(s, dx, out)
            integer, intent(in)                  :: dx
            real(wp), intent(in),   dimension(:) :: s
            real(wp), intent(out),  dimension(:) :: out

            integer                              :: ii

            out = 0.0_wp
            do ii = 1, size(xy)
                out = out + basis_1D(s, Vinv(:, ii), dx) * xy(ii)
            enddo

            return
        end subroutine XorJ
    end function integrate_basis_1d_Ie