using Arpack
include("../ho_basis.jl")

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

Λ = 0
m = 1.0 / 43.281307
V_of_r(r) = -1.2343566 * exp(-r^2 / 100)

# three-body potential V3b = V3b_of_r(r12) * V3b_of_r(r23) * V3b_of_r(r31)
V3b_of_r(r) = 0.43 * exp(-(2/3) * r^2 / 100)

E_max = 30
μω_global = 0.1

μ1ω1 = μω_global * 1/2
basis = ho_basis_2B(E_max, Λ)

V1_elem(l, n1, n2) = V_numerical(V3b_of_r, l, n1, n2; μω_gen=μ1ω1, atol=atol, maxevals=maxevals)
@time "V1" V1 = get_sp_V_matrix(V1_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, E_max=basis.E_max)

V_relative_elem(l, n1, n2) = V_numerical(V3b_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
V_relative_cache = prealloc_V_cache(basis.E_max, ComplexF64)

V_relative_1 = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_relative_cache)
V_relative_2 = get_sp_V_matrix(V_relative_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_relative_cache)
@time "V_relative" V_relative = V_relative_1 * V_relative_2 # commutative?
U = Moshinsky_transform(basis)
@time "V2" V2 =  transpose(U) * V_relative * U

@time "V3b" V3b = V2 * V1 # 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)