# Gaussian potentials in HO space
gl(R, l, n) = throw("unimplemented")

function get_sp_basis(E_max)
    Es = Int[]
    ns = Int[]
    ls = Int[]

    # Heyde p67 with E = N and n = k + 1
    for E in 0 : E_max
        for n in 0 : E ÷ 2
            l = E - 2*n
            push!(Es, E)
            push!(ns, n)
            push!(ls, l)
        end
    end

    return (Es, ns, ls)
end

function get_2p_basis(E_max)
    Es = Int[]
    n1s = Int[]
    l1s = Int[]
    n2s = Int[]
    l2s = Int[]

    # E = 2*n1 + l1 + 2*n2 + l2
    for E in 0 : 2*E_max
        for n1 in 0 : E ÷ 2
            for n2 in 0 : (E - 2*n1) ÷ 2
                for l1 in 0 : (E - 2*n1 - 2*n2)
                    l2 = E - 2*n1 - 2*n2 - l1
                    push!(Es, E)
                    push!(n1s, n1)
                    push!(l1s, l1)
                    push!(n2s, n2)
                    push!(l2s, l2)
                end
            end
        end
    end

    return (Es, n1s, l1s, n2s, l2s)
end

get_V_matrix(V_l, ls, ns) = throw("unimplemented")

get_T_matrix(ns, ls) = throw("unimplemented")

get_H_matrix(V_l, ns, ls) = get_T_matrix(ns, ls) + get_V_matrix(V_l, ns, ls)