using Plots
include("p_space.jl")

# contour
p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])

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

# generating a Berggren basis with a pole using the same system
basis_c = 0.6
basis_E, berg_basis = eigen(get_H_matrix(V_system(basis_c), p, w); permute=false, scale=false)
basis_p = sqrt.(basis_E)
N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis))
berg_basis = berg_basis ./ transpose(N_berg)
berg_basis_w = berg_basis .* w

training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)

training_E = Vector{ComplexF64}(undef, length(training_points))
EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))

# training
for (j, c) in enumerate(training_points)
    H = get_H_matrix(V_system(c), p, w)
    H_berg = transpose(berg_basis_w) * H * berg_basis
    evals, evecs = eigen(H_berg)
    i = argmin(real.(evals))
#    i = identify_pole_i(basis_p, evals)
    training_E[j] = evals[i]
    EC_basis[:, j] = evecs[:, i]
end

EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
N_EC = transpose(EC_basis) * EC_basis

extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)

exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))

# extrapolating
for (j, c) in enumerate(extrapolate_points)
    exact_E[j] = quick_pole_E(V_system(c))

    H = get_H_matrix(V_system(c), p, w)
    H_berg = transpose(berg_basis_w) * H * berg_basis
    H_EC = transpose(EC_basis) * H_berg * EC_basis
    evals = eigvals(H_EC, N_EC)
    i = argmin(abs.(evals .- exact_E[j]))
    extrapolate_E[j] = evals[i]
end

scatter(real.(training_E), imag.(training_E), label="training")
scatter!(real.(exact_E), imag.(exact_E), label="exact")
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
scatter!(real.(basis_E), imag.(basis_E), m=:x, label="Berggren basis")
xlims!(-0.3,0.3)
ylims!(-0.120,0.020)
savefig("temp/" * string(rand(UInt16)) * ".pdf")