From b4e7897cb7532f459d58fd88bceae1a51ed8210e Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Tue, 27 Feb 2024 13:45:43 -0500
Subject: [PATCH] Introduced omega

---
 ho_basis.jl           | 8 ++++----
 ho_basis_benchmark.jl | 5 +++--
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/ho_basis.jl b/ho_basis.jl
index de9f2f3..9fc71e1 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -6,8 +6,8 @@ include("helper.jl")
 inv_factorial(n) = Iterators.prod(inv.(1:n))
 sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
 N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
-Talmi(l, R, k1, k2) = (-1)^(k1 + k2) * (1 + 1/R^2)^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+Talmi(l, R, k1, k2; ω=1.0) = (-1)^(k1 + k2) * (1 + 1/(ω * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
+V_Gaussian(R, l, n1, n2; ω=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; ω=ω) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
 
 function get_sp_basis(E_max)
     Es = Int[]
@@ -55,7 +55,7 @@ end
 
 get_V_matrix(V_l, ls, ns) = throw("unimplemented")
 
-function sp_T_matrix(ns, ls)
+function sp_T_matrix(ns, ls; ω=1.0)
     mat = zeros(length(ns), length(ns))
     for idx in CartesianIndices(mat)
         (i, j) = Tuple(idx)
@@ -70,7 +70,7 @@ function sp_T_matrix(ns, ls)
             end
         end
     end
-    return mat
+    return ω .* mat
 end
 
 get_H_matrix(V_l, ns, ls) = get_T_matrix(ns, ls) + get_V_matrix(V_l, ns, ls)
diff --git a/ho_basis_benchmark.jl b/ho_basis_benchmark.jl
index 42fe176..e43583c 100644
--- a/ho_basis_benchmark.jl
+++ b/ho_basis_benchmark.jl
@@ -2,6 +2,7 @@ using LinearAlgebra, Plots
 include("ho_basis.jl")
 include("p_space.jl")
 
+ω = 0.5
 μ = 0.5
 l = 0
 V1 = -5
@@ -13,8 +14,8 @@ n_max = 15
 ns = collect(0:n_max)
 ls = fill(l, n_max + 1)
 
-T = 1/μ * sp_T_matrix(ns, ls)
-V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns)) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns))
+T = 1/μ * sp_T_matrix(ns, ls; ω=ω)
+V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω)
 
 cs = range(1.35, 0.9, 5)