using SparseArrays, Arpack
include("helper.jl")
include("p_space.jl")

E_target = -0.3919

μ = 0.5
Va = -2
Ra = 2
V_of_r(r) = Va * exp(-r^2 / Ra^2)

vertices = [0, 0.5 - 0.3im, 1, 4]
subdivisions = [16, 16, 16]
ks, ws = get_mesh(vertices, subdivisions)
ls = collect(0:4)

V_l(l, k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)

H_l = [get_H_matrix((k, kp) -> V_l(l, k, kp), ks, ws, μ) for l in ls]

H1 = blockdiag(sparse.(H_l)...)

H = H1
@time "Eigenvalues" evals, _ = eigs(H, nev=10, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
E = nearest(evals, E_target)