using SparseArrays
include("math.jl")

"berg_bases1/2 are lists of (1+l_max) matrices containing the eigenbases corresponding to 1st and 2nd DOFs respectively,
js are a list of tuples (j1, j2) corresponding to 1st and 2nd DOFs respectively,
and ws are the weights needed to evaluate the inner products"
function get_2p_p1p2_matrix(mesh_size, js, Λ, berg_bases1::Vector{Matrix{ComplexF64}}, berg_bases2::Vector{Matrix{ComplexF64}}, ws::Vector{ComplexF64})
    # TODO: Cache / precalculate
    integral1(np, lp, n, l) = sum(berg_bases1[1 + lp][:, np] .* ws .* berg_bases1[1 + l][:, n])
    integral2(np, lp, n, l) = sum(berg_bases2[1 + lp][:, np] .* ws .* berg_bases2[1 + l][:, n])

    basis = iter_prod(js, 1:mesh_size, 1:mesh_size)
    mat = zeros(ComplexF64, length(basis), length(basis))
    Threads.@threads for idx in CartesianIndices(mat)
        (ip, i) = Tuple(idx)
        ((j1p, j2p), n1p, n2p) = basis[ip]
        ((j1, j2), n1, n2) = basis[i]
        val = racahs_reduction_formula(n1p, j1p, n2p, j2p, n1, j1, n2, j2, Λ, integral1, integral2)
        if !(val ≈ 0); mat[idx] = val; end
    end
    return sparse(mat)
end