yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Overflow safing and optimization

This commit is contained in:
Nuwan Yapa 2024-11-19 16:08:16 -05:00
parent 6268f2a7fd
commit 26183767b8
2 changed files with 7 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ho_basis.jl
View File

@ -5,9 +5,10 @@ include("helper.jl")
include("math.jl") include("math.jl")
function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7) 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) const_part = sqrt(μω_gen) * ho_basis_const(l, n1) * ho_basis_const(l, n2)
integrand(r) = ho_basis_func(l, n1, sqrt(μω_gen) * r) * ho_basis_func(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
(integral, _) = quadgk(integrand, 0, Inf; atol=atol, maxevals=maxevals) (integral, _) = quadgk(integrand, 0, Inf; atol=atol, maxevals=maxevals)
return integral return const_part * integral
end end
function get_sp_basis(E_max) function get_sp_basis(E_max)

6
math.jl
View File

@ -12,13 +12,15 @@ V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n
sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1)) sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
sqrt_sqrt_pi = sqrt(sqrt(pi)) 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]) 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) ho_basis_const(l, n) = (-1)^n / sqrt_sqrt_pi * (2.0)^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1)
ho_basis_func(l, n, x) = x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
ho_basis(l, n, x) = ho_basis_const(l, n) * ho_basis_func(l, n, x)
# for implementation of simple relative coordinates # for implementation of simple relative coordinates
double_factorial(n::Int) = Iterators.prod(big, n:-2:1) double_factorial(n::Int) = Iterators.prod(big, n:-2:1)
"Gaussian integral for n ∈ Integers (Ref: Wolfram MathWorld + simplifications)" "Gaussian integral for n ∈ Integers (Ref: Wolfram MathWorld + simplifications)"
gauss_int(a, n) = double_factorial(n - 1) / (2 * a)^((n + 1)/2) * (iseven(n) ? sqrt(π / 2) : 1) gauss_int(a, n) = double_factorial(n - 1) / (2.0 * a)^((n + 1)/2) * (iseven(n) ? sqrt(π / 2) : 1)
"Gives ∫dp p u' u where u' and u' may have different l (Ref: worked out in Mathematica)" "Gives ∫dp p u' u where u' and u' may have different l (Ref: worked out in Mathematica)"
function integral(np, lp, n, l, μω) function integral(np, lp, n, l, μω)