Merge pull request #24 from heltonmc/large_arg_bessely

Output both bessely and besselj for large arguments

Author: heltonmc

src/asymptotics.jl

@@ -1,3 +1,26 @@
+function besseljy_large_argument(v, x::T) where T
+    # gives both (besselj, bessely) for x > 1.6*v
+    α, αp = _α_αp_asymptotic(v, x)
+    b = SQ2OPI(T) / sqrt(αp * x)
+
+    # we need to calculate sin(x - v*pi/2 - pi/4) and cos(x - v*pi/2 - pi/4)
+    # For improved accuracy this is expanded using the formula for sin(x+y+z)
+
+    S, C = sincos(PIO2(T) * v)
+    Sα, Cα = sincos(α)
+
+    CMS = C - S
+    CPS = C + S
+
+    s1 = CMS * Cα
+    s2 = CPS * Sα
+
+    s3 = CMS * Sα
+    s4 = CPS * Cα
+
+    return SQ2O2(T) * (s1 + s2) * b, SQ2O2(T) * (s3 - s4) * b
+end
+
 function _α_αp_asymptotic(v, x::Float64)
     if x > 5*v
         return _α_αp_poly_10(v, x)

src/besselj.jl

@@ -163,7 +163,7 @@ function _besselj(nu, x)
     x < 4.0 && return besselj_small_arguments_orders(nu, x)
 
     large_arg_cutoff = 1.65*nu
-    (x > large_arg_cutoff && x > 20.0) && return besselj_large_argument(nu, x)
+    (x > large_arg_cutoff && x > 20.0) && return besseljy_large_argument(nu, x)[1]
 
 
     debye_cutoff = 2.0 + 1.00035*x + Base.Math._approx_cbrt(302.681*Float64(x))
@@ -190,8 +190,8 @@ function _besselj(nu, x)
     if (debye_diff > large_arg_diff && x > 20.0)
         nu_shift = ceil(Int, large_arg_diff)
         v2 = nu - nu_shift
-        jnu = besselj_large_argument(v2, x)
-        jnum1 = besselj_large_argument(v2 - 1, x)
+        jnu = besseljy_large_argument(v2, x)[1]
+        jnum1 = besseljy_large_argument(v2 - 1, x)[1]
         return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[2]
     else
         nu_shift = ceil(Int, debye_diff)
@@ -203,20 +203,6 @@ function _besselj(nu, x)
     end
 end
 
-# for moderate size arguments of x and v this has relative errors ~9e-15
-# for large arguments relative errors ~1e-13
-function besselj_large_argument(v, x::T) where T
-    α, αp = _α_αp_asymptotic(v, x)
-    b = SQ2OPI(T) / sqrt(αp * x)
-
-    S, C = sincos(PIO2(T)*v)
-    Sα, Cα = sincos(α)
-    s1 = (C - S) * Cα
-    s2 = (C + S) * Sα
-
-    return SQ2O2(T) * (s1 + s2) * b
-end
-
 # generally can only use for x < 4.0
 # this needs a better way to sum these as it produces large errors
 # only valid in non-oscillatory regime (v>1/2, 0<t<sqrt(v^2 - 0.25))

test/bessely_test.jl

@@ -64,3 +64,7 @@ y1_32 = bessely1.(Float32.(x))
 # test that Inf inputs go to zero
 @test bessely1(Inf32) == zero(Float32)
 @test bessely1(Inf64) == zero(Float64)
+
+
+# briefly test the large argument is working
+@test Bessels.besseljy_large_argument(10.0, 100.0)[2] ≈ SpecialFunctions.bessely(10.0, 100.0)