getDoubleMatrixElement Subroutine

    elemental subroutine getDoubleMatrixElement(step1, step2, deltaB, genFlag1, genFlag2, &
                                      bValue, order, x0_element, x1_element)
        ! access the necessary double excitation product terms for the H matrix
        ! element calculation
        !
        ! input:
        ! step1/2   ...     stepvector values if given CSFs, step1 form <d'|
        ! deltaB    ...     current delta b value of excitation
        ! genFlag1/2...     generator types involved in excitation
        ! bValue    ...     current b value of CSF
        ! order     ...     order parameter determining the sign of x1 element
        !
        ! output:
        ! x0_element...     x0 product term (optional)
        ! x1_element...     x1 product term (optional)
        ! -> call this function by x0_element=var, x1_element=var2, if only
        ! specific matrix element is wanted. x0_element is default if only
        ! one output variable is present
        integer, intent(in) :: step1, step2, deltaB, genFlag1, genFlag2
        real(dp), intent(in) :: bValue, order
        real(dp), intent(out), optional :: x0_element, x1_element
        integer :: x0_ind, x1_ind

        if (present(x0_element)) then
            ! first get correct indices to access procedure pointer array
            x0_ind = indArrTwoX0(step1, step2, 3 * deltaB + (genFlag1 + genFlag2) / 2)

            ! then call corresponding function
            x0_element = doubleMatEleX0GUGA(x0_ind)%ptr(bValue)
        end if

        ! same for x1 element
        if (present(x1_element)) then
            x1_ind = indArrTwoX1(step1, step2, 3 * deltaB + (genFlag1 + genFlag2) / 2)

            x1_element = order * doubleMatEleX1GUGA(x1_ind)%ptr(bValue)
        end if

    end subroutine getDoubleMatrixElement