From 2c53ab1f69d962757418ed4b1d03a30ca61c0138 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Wed, 21 Feb 2024 14:40:54 -0500
Subject: [PATCH] Moshinsky matrix (untested)

---
 ho_basis.jl | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/ho_basis.jl b/ho_basis.jl
index 6a6a320..eacedfe 100644
--- a/ho_basis.jl
+++ b/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