add debye approximations

Author: heltonmc

src/U_polynomials.jl

@@ -26,6 +26,20 @@ 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
+    u0 = one(T)
+    u1 = -1 / 24 * evalpoly(p2, (3, -5))
+    u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
+    u3 = -1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425))
+    u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725))
+    u5 = -1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875))
+    u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875))
+    u7 = -1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625,  -221849150488590625))
+    u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125))
+    u9 = -1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375))
+    u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625))
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10))
+end
 
 #=
 function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat

src/besselj.jl

@@ -220,3 +220,16 @@ function log_besselj_small_arguments_orders(v, x::T) where T
     logout = -loggamma(v + 1) + v * log(x/2) + log(out)
     return logout
 end
+
+# valid when x < v (uniform asymptotic expansions)
+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)
+    p = v / sqrt(v^2 - x^2)
+    p2  = v^2/fma(max(v,x), max(v,x), -min(v,x)^2)
+
+    return coef * Uk_poly_Jn(p, v, p2, T)
+end