diff --git a/berggren.jl b/berggren.jl
new file mode 100644
index 0000000..27162c0
--- /dev/null
+++ b/berggren.jl
@@ -0,0 +1,22 @@
+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
diff --git a/berggren_3body_resonance.jl b/berggren_3body_resonance.jl
index 0f7cbf4..57f627b 100644
--- a/berggren_3body_resonance.jl
+++ b/berggren_3body_resonance.jl
@@ -2,6 +2,7 @@ using LinearAlgebra, SparseArrays, Arpack
 include("helper.jl")
 include("p_space.jl")
 include("ho_basis.jl")
+include("berggren.jl")
 
 println("No of threads = ", Threads.nthreads())
 atol = 10^-5
@@ -10,8 +11,7 @@ R_cutoff = 16
 
 Λ = 0
 m = 1.0
-μ1 = m * 1/2
-μ2 = m * 2/3
+μ = m/2 # due to simple relative coordinates
 
 target = 4.0766890719636875 - 0.012758927741074495im
 
@@ -35,7 +35,7 @@ println("Basis size = $basis_size")
 berg_bases = Vector{Matrix{ComplexF64}}(undef, jmax + 1)
 berg_Es = Vector{Vector{ComplexF64}}(undef, jmax + 1)
 for j in 0:jmax
-    berg_E, berg_basis = eigen(get_H_matrix((k, kp) -> V_l(j, k, kp), ks, ws); permute=false, scale=false)
+    berg_E, berg_basis = eigen(get_H_matrix((k, kp) -> V_l(j, k, kp), ks, ws, μ); permute=false, scale=false)
     N_berg = diag(transpose(berg_basis .* ws) * berg_basis)
     berg_basis = berg_basis ./ transpose(sqrt.(N_berg))
     berg_bases[1 + j] = berg_basis
@@ -49,29 +49,31 @@ to_berg_basis(mat, j) = transpose(berg_bases[1 + j] .* ws) * mat * berg_bases[1
 end
 
 @time "T" begin
-    T_blocks = [kron_sum(to_berg_basis(get_T_matrix(ks, μ1), j1), to_berg_basis(get_T_matrix(ks, μ2), j2)) for (j1, j2) in js]
+    T_blocks = [kron_sum(to_berg_basis(get_T_matrix(ks, μ), j1), to_berg_basis(get_T_matrix(ks, μ), j2)) for (j1, j2) in js]
     T = blockdiag(sparse.(T_blocks)...)
 end
 
-@time "V1" begin
-    V1_blocks = [kron(to_berg_basis(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), j1), I(length(ks))) for (j1, _) in js]
-    V1 = blockdiag(sparse.(V1_blocks)...)
+@time "T_cross" T_cross = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ws) ./ (2*μ)
+
+@time "V1 and V2" begin
+    V_blocks = [kron_sum(to_berg_basis(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), j1), to_berg_basis(get_V_matrix((k, kp) -> V_l(j2, k, kp), ks, ws), j2)) for (j1, j2) in js]
+    V = blockdiag(sparse.(V_blocks)...)
 end
 
 E_max = 30
 μω_global = 0.5
-μ1ω1 = μω_global * 1/2
-μ2ω2 = μω_global * 2
+# due to simple relative coordinates
+μω = μω_global * 2
+μ = m/2
 
-@time "V2_HO" V2_HO =  get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global)
+@time "V12_HO" V12_HO =  get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
 
-@time "W_right" W_right = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=true)
-@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=true)
+@time "W" W = get_W_matrix(basis, E_max, Λ, μω, μω; weights=true)
 
-@time "V2_p" V2_p = W_left * V2_HO * transpose(W_right)
-@time "V2" V2 = transpose(U) * V2_p * U
+@time "V12_p" V12_p = W * V12_HO * transpose(W)
+@time "V12" V12 = transpose(U) * V12_p * U
 
-@time "H" H = T + V1 + V2
+@time "H" H = T + T_cross + V + V12
 
 @time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
 display(evals)