using LinearAlgebra, Plots
include("ho_basis.jl")
include("p_space.jl")

μω_gen = 0.5 * exp(-1im * 0.45 * pi)
μ = 0.5
l = 0
V1 = -5
R1 = sqrt(3)
V2 = 2
R2 = sqrt(10)
n_max = 15

ns = collect(0:n_max)
ls = fill(l, n_max + 1)

T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)

V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; μω_gen=μω_gen) + V2 * V_Gaussian(R2, l, n1, n2; μω_gen=μω_gen)
V = get_sp_V_matrix(V_l, ns, ls; dtype=ComplexF64)

cs = range(1.25, 0.25, 10)

E = similar(cs, ComplexF64)
bench_E = similar(cs, ComplexF64)

for (j, c) in enumerate(cs)
    H = T + c .* V
    evals = eigvals(collect(H))
    bench_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.4, meshpoints=512)
    i = argmin(abs.(evals .- bench_E[j]))
    E[j] = evals[i]
end

scatter(real.(bench_E), imag.(bench_E), label="p-space")
scatter!(real.(E), imag.(E), label="HO basis")
xlims!(-0.7,0.4)
ylims!(-0.25,0.25)