From 98e8cdb9d3fbfc164a1c9b4d24c5169ac1b7c207 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Fri, 26 Apr 2024 17:02:32 -0400
Subject: [PATCH] Numerical V matrix elements

---
 p_space.jl | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/p_space.jl b/p_space.jl
index e63c854..7671c43 100644
--- a/p_space.jl
+++ b/p_space.jl
@@ -1,9 +1,5 @@
 using LinearAlgebra
-using FastGaussQuadrature
-
-# Gaussian potentials in momentum space
-g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))
-g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)
+using SpecialFunctions, FastGaussQuadrature, QuadGK
 
 function gausslegendre_shifted(a, b, n)
     scale = (b - a) / 2
@@ -50,4 +46,16 @@ function quick_pole_E(V_pq, μ=0.5; cs_angle=0.4, cutoff=8.0, meshpoints=256)
     p, w = get_mesh([0, cutoff * exp(-1im * cs_angle)], [meshpoints])
     evals = eigvals(get_H_matrix(V_pq, p, w, μ))
     return evals[identify_pole_i(p, evals, μ)]
-end
\ No newline at end of file
+end
+
+# Gaussian potentials in momentum space
+g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))
+g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)
+
+# general potential (numerical integration)
+jHat(l, z) = z * sphericalbesselj(l, z)
+function Vl_mat_elem(V_of_r, l, p, q; atol=0, maxevals=10^7, R_cutoff=Inf)
+    integrand(r) = jHat(l, p * r) * V_of_r(r) * jHat(l, q * r)
+    (integral, _) = quadgk(integrand, 0, R_cutoff; atol=atol, maxevals=maxevals)
+    return (2 / pi) * integral
+end