Merge branch 'ho_basis' into ho_basis_log

ho_basis.jl

@@ -3,11 +3,9 @@ using SpecialFunctions
 include("helper.jl")
 
 # Gaussian potentials in HO space
-inv_factorial(n) = Iterators.prod(inv.(1:n))
+log_N(l, n, k) = -logfactorial(n)/2 + logabsbinomial(n, k)[1] + loggamma(n + l + 3/2)/2 - loggamma(k + l + 3/2)
-sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
+log_Talmi(l, R, k1, k2) = -(3/2 + l + k1 + k2) * log(1 + 1/R^2) + loggamma(3/2 + l + k1 + k2)
-N_lnk(l, n, k) = (-1)^(n-k) * inv_factorial(k) * binomial(n+l+1/2, n-k) * sqrt_factorial(n) / sqrt(gamma(n+l+3/2))
+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)])
-Talmi(l, R, k1, k2) = (1 + 1/R^2)^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2) = better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
 
 function get_sp_basis(E_max)
     Es = Int[]