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

berggren_mesh = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
csm_mesh = get_mesh([0, 8 - 3im], [512])

for mesh in (berggren_mesh, csm_mesh)
    p, w = mesh
    mesh_E = p.*p ./ (2*0.5)

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

    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))

    for (j, c) in enumerate(training_points)
        evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
        i = identify_pole_i(p, evals)
        training_E[j] = evals[i]
        EC_basis[:, j] = evecs[:, i]
    end

    EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
    EC_basis_w = EC_basis .* w
    N_EC = transpose(EC_basis_w) * 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))

    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_EC = transpose(EC_basis_w) * H * 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")
    plot!(real.(mesh_E), imag.(mesh_E), label="contour")
    xlims!(-0.3,0.3)
    ylims!(-0.120,0.020)
    savefig("temp/" * string(rand(UInt16)) * ".pdf")
end