diff --git a/ho_basis.jl b/ho_basis.jl
index d6c9a8c..63a0d3a 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -5,9 +5,10 @@ include("helper.jl")
 include("math.jl")
 
 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)
-    return integral
+    return const_part * integral
 end
 
 function get_sp_basis(E_max)
diff --git a/math.jl b/math.jl
index 03b1869..d407979 100644
--- a/math.jl
+++ b/math.jl
@@ -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_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)
+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
 double_factorial(n::Int) = Iterators.prod(big, n:-2:1)
 
 "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)"
 function integral(np, lp, n, l, μω)