23 lines
1.1 KiB
Julia
23 lines
1.1 KiB
Julia
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
|