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

Merge branch 'ho_basis' into ho_basis_log

This commit is contained in:
Nuwan Yapa 2024-02-27 11:44:46 -05:00
parent 0d9d5bd44d 5d39015976
commit 8e8880f427
2 changed files with 12 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
helper.jl Normal file
View File

@ -0,0 +1,10 @@
"Sum over array while minimizing catastrophic cancellation as much as possible"
function better_sum(arr::Array{Float64})
pos_arr = arr[arr .> 0]
neg_arr = arr[arr .< 0]
sort!(pos_arr)
sort!(neg_arr, rev=true)
return sum(pos_arr) + sum(neg_arr)
end

3
ho_basis.jl
View File

@ -1,10 +1,11 @@
using NuclearToolkit
using SpecialFunctions
include("helper.jl")
# Gaussian potentials in HO space
log_N(l, n, k) = -logfactorial(n)/2 + logabsbinomial(n, k)[1] + loggamma(n + l + 3/2)/2 - loggamma(k + l + 3/2)
log_Talmi(l, R, k1, k2) = -(3/2 + l + k1 + k2) * log(1 + 1/R^2) + loggamma(3/2 + l + k1 + k2)
V_Gaussian(R, l, n1, n2) = (-1)^(n1 + n2) * sum([(-1)^(k1 + k2) * exp(log_N(l, n1, k1) + log_N(l, n2, k2) + log_Talmi(l, R, k1, k2)) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
V_Gaussian(R, l, n1, n2) = (-1)^(n1 + n2) * better_sum([(-1)^(k1 + k2) * exp(log_N(l, n1, k1) + log_N(l, n2, k2) + log_Talmi(l, R, k1, k2)) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
function get_sp_basis(E_max)
Es = Int[]