using Plots
include("../helper.jl")
include("../p_space.jl")

vertices = [0, 0.5 - 0.3im, 1, 6]
subdivisions = [64, 64, 64]
p, w = get_mesh(vertices, subdivisions)

# resonance example from my thesis
V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)

basis_eig = eigen(get_H_matrix(V_basis, p, w))
basis = basis_eig.vectors .* w

basis_evals = basis_eig.values
E_target = 0.7
E = nearest(basis_evals, E_target)

print("pole E = $E")

# ResonanceEC: Eq. (20)
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))

H = get_H_matrix(V_system(0.45), p, w)
H_berggren = transpose(basis) * H * basis

evals = eigvals(H)
scatter(real.(evals), imag.(evals), label="E")
scatter!(real.(basis_evals), imag.(basis_evals), label="basis", markershape=:x)
xlims!((0, 1))