From 71c789376dcde6c7767d50f67a98f279d0a17c50 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Mon, 22 Apr 2024 15:12:16 -0400
Subject: [PATCH] Numerical calculation of V matrix

---
 ho_basis.jl         | 13 +++++++++++++
 test/numerical_V.jl | 20 ++++++++++++++++++++
 2 files changed, 33 insertions(+)
 create mode 100644 test/numerical_V.jl

diff --git a/ho_basis.jl b/ho_basis.jl
index 268e448..b97576e 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,5 +1,6 @@
 using SparseArrays
 using SpecialFunctions
+using QuadGK
 include("helper.jl")
 
 # Gaussian potentials in HO space
@@ -120,3 +121,15 @@ function pick_Moshinsky_bracket(BRAC, n1′, l1′, n2′, l2′, n1, l1, n2, l2
     # BRAC(NP,N1P,MP,N1,N2,N,M,L)
     return BRAC[1 + NP, 1 + n1′, 1 + MP, 1 + n1, 1 + n2, 1 + N, 1 + M, 1]
 end
+
+# numerical evaluation of V matrix elements
+sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
+sqrt_sqrt_pi = sqrt(sqrt(pi))
+laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k)  for k in 0:n])
+ho_basis(l, n, x) = (-1)^n / sqrt_sqrt_pi * 2^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1) * x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
+
+function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0)
+    integrand(r) = sqrt(μω_gen) * ho_basis(l, n1, sqrt(μω_gen) * r) * ho_basis(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
+    (integral, _) = quadgk(integrand, 0, Inf)
+    return integral
+end
diff --git a/test/numerical_V.jl b/test/numerical_V.jl
new file mode 100644
index 0000000..b419e45
--- /dev/null
+++ b/test/numerical_V.jl
@@ -0,0 +1,20 @@
+using Statistics
+
+include("../ho_basis.jl")
+
+ls = 0:4
+ns = 0:15
+
+R = 2
+μω_gen = 0.3
+V_of_r(r) = exp(-r^2 / R^2)
+
+for l in ls
+    println("Testing l=", l)
+
+    @time "Analytical" V_ana = [V_Gaussian(R, l, n1, n2; μω_gen=μω_gen) for (n1, n2) in Iterators.product(ns, ns)]
+    @time "Numerical" V_num = [V_numerical(V_of_r, l, n1, n2; μω_gen=μω_gen) for (n1, n2) in Iterators.product(ns, ns)]
+
+    mean_diff = Statistics.mean(abs.((V_num .- V_ana) ./ V_ana))
+    @assert mean_diff < 1e-5 "Mean absoute difference of $(mean_diff * 100)% between analytical and numerical calculations"
+end