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

function get_sp_basis(max_N)
    ls = Int[]
    ks = Int[]

    # Heyde p67
    for N in max_N:-1:0
        for l in N:-2:0
            k = (N - l) / 2
            push!(ls, l)
            push!(ks, k)
        end
    end

    return (ls, ks)
end

get_V_matrix(V_lk, ls, ks) = throw("unimplemented")

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

get_H_matrix(V_ls, ls, ks) = get_T_matrix(ls, ks) + get_V_matrix(V_lk, ls, ks)