update coefficients

Author: heltonmc

src/U_polynomials.jl

@@ -26,7 +26,7 @@ function Uk_poly_In(p, v, p2, ::Type{Float32})
     u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
     return evalpoly(-p/v, (u0, u1, u2))
  end
- function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
+ function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Union{Float64, BigFloat}
     u0 = one(T)
     u1 = -1 / 24 * evalpoly(p2, (3, -5))
     u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))

src/besselj.jl

@@ -226,10 +226,11 @@ function besselj_debye(v, x)
     T = eltype(x)
     S = promote_type(T, Float64)
     x = S(x)
-    n = v*log(v/x + sqrt((v/x)^2 - 1)) - sqrt(v^2 - x^2)
-    coef = inv(sqrt(2*T(pi))) * exp(-n) / (v^2 - x^2)^(1/4)
+    b = sqrt(v^2 - x^2)
+    n = v * log(v/x + sqrt((v/x)^2 - 1)) - b
+    coef = inv(sqrt(2*T(pi))) * exp(-n) / sqrt(b)
     p = v / sqrt(v^2 - x^2)
-    p2  = v^2/fma(max(v,x), max(v,x), -min(v,x)^2)
+    p2  = v^2/fma(v,v, -x^2)
 
     return coef * Uk_poly_Jn(p, v, p2, T)
 end