From 6e5e1d013282221d15a6de4c25c4ad48a7fa4a60 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Thu, 19 Sep 2024 18:16:37 -0400 Subject: [PATCH] Reorganize --- ho_basis.jl | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/ho_basis.jl b/ho_basis.jl index b7826b3..0d1bcf0 100644 --- a/ho_basis.jl +++ b/ho_basis.jl @@ -11,6 +11,20 @@ N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2) V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)]) +# 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, atol=0, maxevals=10^7) + 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; atol=atol, maxevals=maxevals) + return integral +end + +############################################################## + function get_sp_basis(E_max) Es = Int[] ns = Int[] @@ -126,18 +140,6 @@ function pick_Moshinsky_bracket(BRAC, n1′, l1′, n2′, l2′, n1, l1, n2, l2 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, atol=0, maxevals=10^7) - 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; atol=atol, maxevals=maxevals) - return integral -end - function get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global; atol=10^-6, maxevals=10^5) _, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ) l_max = max(maximum(l1s), maximum(l2s))