Moshinsky matrix (untested)

ho_basis.jl

@@ -1,3 +1,5 @@
+using NuclearToolkit
+
 # Gaussian potentials in HO space
 gl(R, l, n) = throw("unimplemented")
 
@@ -50,3 +52,27 @@ get_V_matrix(V_l, ls, ns) = throw("unimplemented")
 get_T_matrix(ns, ls) = throw("unimplemented")
 
 get_H_matrix(V_l, ns, ls) = get_T_matrix(ns, ls) + get_V_matrix(V_l, ns, ls)
+
+function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
+    l_max = 2*max(maximum(l1s), maximum(l2s)) # OPTIMIZE
+    E_max = maximum(Es) # OPTIMIZE
+    j_max = l_max # OPTIMIZE
+    to = 0 # unused
+    
+    dtri = NuclearToolkit.prep_dtri(l_max) ;     println("dtri prepared")
+    dcgm0 = NuclearToolkit.prep_dcgm0(l_max) ;     println("dcgm0 prepared")
+    d6j = NuclearToolkit.prep_d6j_int(E_max, j_max, to) ;     println("d6j prepared")
+
+    mat = zeros(length(Es), length(Es))
+    s = hcat(Es, n1s, l1s, n2s, l2s)
+    for idx in CartesianIndices(mat)
+        (i, j) = Tuple(idx)
+        (Elhs, N, L, n, l) = s[i, :]
+        (Erhs, n1, l1, n2, l2) = s[j, :]
+        if Elhs == Erhs
+            mat[i, j] = NuclearToolkit.HObracket_d6j(N, L, n, l, n1, l1, n2, l2, Λ, 1.0, dtri, dcgm0, d6j, to)
+        end
+    end
+
+    return mat
+end