using QuadGK
include("../math.jl")

# Testing ∫dp p u' u where u' and u' may have different l

function explicit_integral(np, lp, n, l, μω)
    integrand(p) = p * (-1)^(n + np) * 1/sqrt(μω) * ho_basis(l, n, 1/sqrt(μω) * p) * ho_basis(lp, np, 1/sqrt(μω) * p)
    (integral, _) = quadgk(integrand, 0, Inf; atol=1e-10)
    return integral
end

n_list = 0 : 5
l_list = 0 : 5
μω_list = [1, 2.66]

for np in n_list
    for n in n_list
        for lp in l_list
            for l in l_list
                for μω in μω_list
                    println("Checking n=$n, np=$np, l=$l, lp=$lp, μω=$μω")
                    expl = explicit_integral(np, lp, n, l, μω)
                    anal = integral(np, lp, n, l, μω)
                    @assert isapprox(expl, anal; atol=1e-8) "Explicit=$(expl) and analytical=$(anal)"
                end
            end
        end
    end
end

# Testing the full Racah's reduction formula

test_data = DelimitedFiles.readdlm("test/racah_test_data.csv", ',') # format = n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2, racah_matrix_element
test_data_ints = Int.(test_data[:, 1:(end - 3)])
test_data_floats = test_data[:, (end - 2):end]

for i in 1 : size(test_data)[1]
    n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ = test_data_ints[i, :]
    μ1ω1, μ2ω2, racah_matrix_element = test_data_floats[i, :]
    println("Checking n1'=$n1p, l1'=$l1p, n2'=$n2p, l2'=$l2p, n1=$n1, l1=$l1, n2=$n2, l2=$l2, Λ=$Λ, μ1ω1=$μ1ω1, μ2ω2=$μ2ω2")
    val = racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2)
    @assert isapprox(val, racah_matrix_element; atol=1e-8) "Formula=$(val) and test reference=$racah_matrix_element"
end