From 58e32f4950c950290f32ca61605654207713f62f Mon Sep 17 00:00:00 2001
From: ysyapa <ysyapa@ncsu.edu>
Date: Fri, 31 May 2024 20:13:20 +0000
Subject: [PATCH] Big optimization due to triangle inequality

---
 3body_resonance.jl | 2 +-
 ho_basis.jl        | 3 ++-
 ho_basis_3body.jl  | 2 +-
 test/moshinsky.jl  | 2 +-
 4 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/3body_resonance.jl b/3body_resonance.jl
index 8d8e9ff..97453ec 100644
--- a/3body_resonance.jl
+++ b/3body_resonance.jl
@@ -17,7 +17,7 @@ E_max = 30
 println("No of threads = ", Threads.nthreads())
 
 @time "Basis" begin
-    Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
+    Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
     l_max = max(maximum(l1s), maximum(l2s))
     n_max = max(maximum(n1s), maximum(n2s))
     mask1 = (n2s .== n2s') .&& (l2s .== l2s')
diff --git a/ho_basis.jl b/ho_basis.jl
index 84dda6f..0ec0bb4 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -29,7 +29,7 @@ function get_sp_basis(E_max)
     return (Es, ns, ls)
 end
 
-function get_2p_basis(E_max)
+function get_2p_basis(E_max, Λ=-1)
     Es = Int[]
     n1s = Int[]
     l1s = Int[]
@@ -42,6 +42,7 @@ function get_2p_basis(E_max)
             for n2 in 0 : (E - 2*n1) ÷ 2
                 for l1 in 0 : (E - 2*n1 - 2*n2)
                     l2 = E - 2*n1 - 2*n2 - l1
+                    if Λ≥0 && !triangle_ineq(l1, l2, Λ); continue; end
                     push!(Es, E)
                     push!(n1s, n1)
                     push!(l1s, l1)
diff --git a/ho_basis_3body.jl b/ho_basis_3body.jl
index 57bb6d6..427cdcd 100644
--- a/ho_basis_3body.jl
+++ b/ho_basis_3body.jl
@@ -18,7 +18,7 @@ E_max = 40
 println("No of threads = ", Threads.nthreads())
 
 @time "Basis" begin
-    Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
+    Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
     mask1 = (n2s .== n2s') .&& (l2s .== l2s')
     mask2 = (n1s .== n1s') .&& (l1s .== l1s')
 end
diff --git a/test/moshinsky.jl b/test/moshinsky.jl
index 772b6cd..515a108 100644
--- a/test/moshinsky.jl
+++ b/test/moshinsky.jl
@@ -9,7 +9,7 @@ E_max = 30
 
 println("No of threads = ", Threads.nthreads())
 
-@time "Basis" Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
+@time "Basis" Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
 
 println("Basis size = ", length(Es))