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

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

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

    return (ns, ls)
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)