From 6b7b7599796fde914654679dd39da2407c394c01 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Wed, 21 Feb 2024 19:32:11 -0500
Subject: [PATCH] Gaussian potential implemented

---
 ho_basis.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/ho_basis.jl b/ho_basis.jl
index d97f53a..b79d646 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,7 +1,11 @@
 using NuclearToolkit
+using SpecialFunctions
 
 # Gaussian potentials in HO space
-gl(R, l, n) = throw("unimplemented")
+N_nl(n, l) = sqrt(2*factorial(n) / gamma(l + n + 3/2)^3)
+prefactor(n, l, k) = (-1)^k * binomial(n + l + 1/2, n - k) / factorial(k)
+Talmi(l, R, k1, k2) = (1/2) / (1 + 1/R^2)^(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
+V_Gaussian(R, l, n1, n2) = N_nl(n1, l) * N_nl(n2, l) * sum([prefactor(n1, l, k1) * prefactor(n2, l, k2) * Talmi(l, R, k1, k2) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
 
 function get_sp_basis(E_max)
     Es = Int[]