From 3200b7e0415ae789a82d178979a21759cf4bd16c Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Sat, 4 May 2024 21:40:43 -0400
Subject: [PATCH] EC for 3 body resonance

---
 3body_resonance_EC.jl | 53 +++++++++++++++++++++++++++++++++++++++++++
 helper.jl             |  5 +++-
 2 files changed, 57 insertions(+), 1 deletion(-)
 create mode 100644 3body_resonance_EC.jl

diff --git a/3body_resonance_EC.jl b/3body_resonance_EC.jl
new file mode 100644
index 0000000..426c18e
--- /dev/null
+++ b/3body_resonance_EC.jl
@@ -0,0 +1,53 @@
+using Plots
+include("3body_resonance.jl")
+
+current_E = 5.9673 - 0.0006im
+
+training_c = 2.0 : -0.2 : 1.2
+extrapolating_c = 1.05 .- [0.0 : 0.1 : 0.4; 0.45 : 0.05 : 0.60]
+
+@time "H0" H0 = T1 + T2
+@time "V" V = V1 + V2
+
+# free memory
+Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
+GC.gc()
+
+exact = ComplexF64[]
+training = ComplexF64[]
+extrapolated = ComplexF64[]
+training_vecs = Vector{ComplexF64}[]
+
+for c in training_c
+    println("Training for c = $c")
+    H = H0 + c .* V
+    evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
+
+    global current_E = nearest(evals, current_E)
+    push!(training, current_E)
+    push!(training_vecs, evecs[:, nearestIndex(evals, current_E)])
+end
+
+EC_basis = hcat(training_vecs...)
+N_EC = transpose(EC_basis) * EC_basis
+H0_EC = transpose(EC_basis) * H0 * EC_basis
+V_EC = transpose(EC_basis) * V * EC_basis
+
+for c in extrapolating_c
+    println("Extrapolating for c = $c")
+    H = H0 + c .* V
+    evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
+
+    global current_E = nearest(evals, current_E)
+    push!(exact, current_E)
+
+    # extrapolation
+    H_EC = H0_EC + c .* V_EC
+    evals = eigvals(H_EC, N_EC)
+    push!(extrapolated, nearest(evals, current_E))
+end
+
+scatter(real.(training),imag.(training), label="training")
+scatter!(real.(exact),imag.(exact), label="exact")
+scatter!(real.(extrapolated),imag.(extrapolated), label="extrapolated")
+savefig("temp/NCSM.pdf")
\ No newline at end of file
diff --git a/helper.jl b/helper.jl
index 96e42e3..0307729 100644
--- a/helper.jl
+++ b/helper.jl
@@ -16,8 +16,11 @@ better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(i
 "The triangle inequality for angular momenta"
 triangle_ineq(l1, l2, L) = abs(l1 - l2) ≤ L && L ≤ (l1 + l2)
 
+"Index of the nearest value in a list to a given reference point"
+nearestIndex(list::Array, ref) = argmin(norm.(list .- ref))
+
 "Nearest value in a list to a given reference point"
-nearest(list::Array, ref) = list[argmin(norm.(list .- ref))]
+nearest(list::Array, ref) = list[nearestIndex(list, ref)]
 
 "Simple implementation of the Kronecker sum"
 function kron_sum(A::AbstractMatrix, B::AbstractMatrix)