using Arpack
include("ho_basis.jl")

coord_system = src
atol = 10 ^ -5
maxevals = 10 ^ 5

Λ = 0
m = 1.0
V_of_r(r) = -2 * exp(-r^2 / 4)

# three-body potential V3b = V3b_of_r(r12) * V3b_of_r(r23) * V3b_of_r(r31)
V3b_of_r(r) = -1 * exp(-r^2 / 4)

E_max = 40
μω_global = 0.3

# compute three-body potential matrix elements
μω = μω_global * 2
basis = ho_basis_2B(E_max, Λ)
V3b_elem(l, n1, n2) = V_numerical(V3b_of_r, l, n1, n2; μω_gen=μω, atol=atol, maxevals=maxevals)
V_cache = prealloc_V_cache(basis.E_max, ComplexF64)
@time "V1" V1 = get_sp_V_matrix(V3b_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_cache)
@time "V2" V2 = get_sp_V_matrix(V3b_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_cache)
@time "V12" V12 = get_src_V12_matrix(V3b_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
@time "V3b" V3b = (V1 * V2) * V12 # commutative?

H = get_3b_H_matrix(coord_system, V_of_r, μω_global, E_max, Λ, m) + V3b

@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
display(evals)