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

Move math functions to a seperate file

This commit is contained in:
Nuwan Yapa 2024-09-19 18:25:30 -04:00
parent 6e5e1d0132
commit c2f7bde00e
2 changed files with 15 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
ho_basis.jl
View File

@ -1,21 +1,8 @@
using SparseArrays using SparseArrays
using SpecialFunctions
using QuadGK using QuadGK
using LRUCache using LRUCache
include("helper.jl") include("helper.jl")
include("math.jl")
# Gaussian potentials in HO space
inv_factorial(n) = Iterators.prod(inv.(1:n))
sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + 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) 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) integrand(r) = sqrt(μω_gen) * ho_basis(l, n1, sqrt(μω_gen) * r) * ho_basis(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
@ -23,8 +10,6 @@ function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
return integral return integral
end end
##############################################################
function get_sp_basis(E_max) function get_sp_basis(E_max)
Es = Int[] Es = Int[]
ns = Int[] ns = Int[]

14
math.jl Normal file
View File

@ -0,0 +1,14 @@
using SpecialFunctions
# Gaussian potentials in HO space
inv_factorial(n) = Iterators.prod(inv.(1:n))
sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + 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)])
# for 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)