add support to compute both Jv and Yv and Hankel

Author: heltonmc

src/Bessels.jl

@@ -42,6 +42,6 @@ include("Polynomials/besselj_polys.jl")
 include("asymptotics.jl")
 include("gamma.jl")
 
-#include("hankel.jl")
+include("hankel.jl")
 
 end

src/U_polynomials.jl

@@ -41,7 +41,7 @@ function besseljy_debye(v, x)
     return coef_Jn * Uk_Jn, coef_Yn * Uk_Yn
 end
 
-besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && x > 15
+besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && nu > 15
 
 """
     hankel_debye(nu, x::T)

src/bessely.jl

@@ -310,10 +310,10 @@ function bessely_positive_args(nu, x::T) where T
     hankel_debye_cutoff(nu, x) && return imag(hankel_debye(nu, x))
 
     # use power series for small x and for when nu > x
-    bessely_series_cutoff(nu, x) && return bessely_power_series(nu, x)
+    bessely_series_cutoff(nu, x) && return bessely_power_series(nu, x)[1]
 
     # for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
-    besseljy_chebyshev_cutoff(nu, x) && return bessely_chebyshev(nu, x)
+    besseljy_chebyshev_cutoff(nu, x) && return bessely_chebyshev(nu, x)[1]
 
     # at this point x > 19.0 (for Float64) and fairly close to nu
     # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
@@ -340,6 +340,7 @@ end
 
 Computes ``Y_{nu}(x)`` using the power series when nu is not an integer.
 In general, this is most accurate for small arguments and when nu > x.
+Outpus both (Y_{nu}(x), J_{nu}(x)).
 """
 function bessely_power_series(v, x::T) where T
     MaxIter = 3000
@@ -358,7 +359,7 @@ function bessely_power_series(v, x::T) where T
         b *= -inv((-v + i + one(T)) * (i + one(T))) * t2
     end
     s, c = sincospi(v)
-    return (out*c - out2) / s
+    return (out*c - out2) / s, out
 end
 bessely_series_cutoff(v, x) = (x < 7.0) || v > 1.35*x - 4.5
 
@@ -375,7 +376,7 @@ Forward recurrence is used to fill orders starting at low orders ν ∈ (0, 2).
 function bessely_chebyshev(v, x)
     v_floor, _ = modf(v)
     Y0, Y1 = bessely_chebyshev_low_orders(v_floor, x)
-    return besselj_up_recurrence(x, Y1, Y0, v_floor + 1, v)[1]
+    return besselj_up_recurrence(x, Y1, Y0, v_floor + 1, v)
 end
 
 # only implemented for Float64 so far

src/hankel.jl

@@ -1,4 +1,134 @@
-#= Aren't tested yet
+function besseljy(nu::Real, x::T) where T
+    isinteger(nu) && return besseljy(Int(nu), x)
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+
+    Ynu = bessely_positive_args(abs_nu, abs_x)
+    Jnu = bessely_positive_args(abs_nu, abs_x)
+
+    if nu >= zero(T)
+        if x >= zero(T)
+            return Jnu, Ynu
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cospi(abs_nu)
+        end
+    else
+        spi, cpi = sincospi(abs_nu)
+        out = Jnu * cpi - Ynu * spi
+        if x >= zero(T)
+            return out, Ynu * cpi + Jnu * spi
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return cpi * (Ynu * cispi(nu) + 2im * Jnu * cpi) + Jnu * spi * cispi(abs_nu)
+        end
+    end
+end
+
+function besseljy(nu::Integer, x::T) where T
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+    sg = iseven(abs_nu) ? 1 : -1
+
+    Jnu = besselj_positive_args(abs_nu, abs_x)
+    Ynu = bessely_positive_args(abs_nu, abs_x)
+    if nu >= zero(T)
+        if x >= zero(T)
+            return Jnu, Ynu
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Jnu * sg, Ynu * sg + 2im * sg * besselj_positive_args(abs_nu, abs_x)
+        end
+    else
+        if x >= zero(T)
+            return Jnu * sg, Ynu * sg
+        else
+            spi, cpi = sincospi(abs_nu)
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return (cpi*Jnu - spi*Ynu) * sg, Ynu + 2im * besselj_positive_args(abs_nu, abs_x)
+        end
+    end
+end
+
+function besseljy_positive_args(nu::Real, x::T) where T
+    nu == 0 && return (besselj0(x), bessely0(x))
+    nu == 1 && return (besselj1(x), bessely1(x))
+    iszero(x) && return (besselj_positive_args(nu, x), bessely_positive_args(nu, x))
+
+    # x < ~nu branch see src/U_polynomials.jl
+    besseljy_debye_cutoff(nu, x) && return besseljy_debye(nu, x)
+
+    # large argument branch see src/asymptotics.jl
+    besseljy_large_argument_cutoff(nu, x) && return besseljy_large_argument(nu, x)
+
+    # x > ~nu branch see src/U_polynomials.jl on computing Hankel function
+    if hankel_debye_cutoff(nu, x)
+        H = hankel_debye(nu, x)
+        return real(H), imag(H)
+    end
+
+    # use forward recurrence if nu is an integer up until the continued fraction becomes inefficient
+    if isinteger(nu) && nu < 150
+        Y0 = bessely0(x)
+        Y1 = bessely1(x)
+        Ynm1, Yn = besselj_up_recurrence(x, bessely1(x), bessely0(x), 1, nu-1)
+
+        ratio_Jv_Jvm1 = besselj_ratio_jnu_jnum1(nu, x)
+        Jn = 2 / (π*x * (Ynm1 - Yn / ratio_Jv_Jvm1))
+        return Jn, Yn
+    end
+
+    # use power series for small x and for when nu > x
+    if bessely_series_cutoff(nu, x) && besselj_series_cutoff(nu, x)
+        Yn, Jn = bessely_power_series(nu, x)
+        return Jn, Yn
+    end
+
+    # for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
+    if besseljy_chebyshev_cutoff(nu, x)
+        Yn, Ynp1 = bessely_chebyshev(nu, x)
+        ratio_Jvp1_Jv = besselj_ratio_jnu_jnum1(nu+1, x)
+        Jnp1 = 2 / (π*x * (Yn - Ynp1 / ratio_Jvp1_Jv))
+        return Jnp1 / ratio_Jvp1_Jv, Yn
+    end
+
+    # at this point x > 19.0 (for Float64) and fairly close to nu
+    # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
+    nu_shift = floor(nu) - ceil(Int, -1.5 + x + Base.Math._approx_cbrt(-411*x))
+    v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
+
+    Hnu = hankel_debye(v2, x)
+    Hnum1 = hankel_debye(v2 - 1, x)
+    # forward recurrence is stable for Hankel when x >= nu
+    if x >= nu
+        H = besselj_up_recurrence(x, Hnu, Hnum1, v2, nu)[1]
+        return real(H), imag(H)
+    else
+        Yn, Ynp1 = besselj_up_recurrence(x, imag(Hnu), imag(Hnum1), v2, nu)
+        ratio_Jvp1_Jv = besselj_ratio_jnu_jnum1(nu+1, x)
+        Jnp1 = 2 / (π*x * (Yn - Ynp1 / ratio_Jvp1_Jv))
+        return Jnp1 / ratio_Jvp1_Jv, Yn
+    end
+end
+
+function besselj_ratio_jnu_jnum1(n, x::T) where T
+    MaxIter = 5000
+    xinv = inv(x)
+    xinv2 = 2 * xinv
+    d = x / (n + n)
+    a = d
+    h = a
+    b = muladd(2, n, 2) * xinv
+    for _ in 1:MaxIter
+        d = inv(b - d)
+        a *= muladd(b, d, -1)
+        h = h + a
+        b = b + xinv2
+        abs(a / h) <= eps(T) && break
+    end
+    return h
+end
+
 function besselh(nu::Float64, k::Integer, x::AbstractFloat)
     if k == 1
         return complex(besselj(nu, x), bessely(nu, x))
@@ -11,4 +141,3 @@ end
 
 hankelh1(nu, z) = besselh(nu, 1, z)
 hankelh2(nu, z) = besselh(nu, 2, z)
-=#

test/besselj_test.jl

@@ -85,6 +85,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besselj(BigFloat(v), BigFloat(xx))
     @test isapprox(besselj(v, xx), Float64(sf), rtol=5e-14)
+    @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-14)
 end
 
 # test half orders (SpecialFunctions does not give big float precision)
@@ -94,15 +95,19 @@ x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97,
 nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235, 51.23, 72.23435, 80.5, 98.5, 104.2]
 for v in nu, xx in x
     xx *= v
-    @test isapprox(besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=1e-12)
+    sf = SpecialFunctions.besselj(v, xx)
+    @test isapprox(besselj(v, xx), sf, rtol=1e-12)
+    @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], sf, rtol=1e-12)
 end
 
 ## test large orders
 nu = [150, 165.2, 200.0, 300.0, 500.0, 1000.0, 5000.2, 10000.0, 50000.0]
 x = [0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92,0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99,0.995, 0.999, 1.0, 1.01, 1.05, 1.08, 1.1, 1.2]
 for v in nu, xx in x
     xx *= v
-    @test isapprox(besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=5e-11)
+    sf = SpecialFunctions.besselj(v, xx)
+    @test isapprox(besselj(v, xx), sf, rtol=5e-11)
+    @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-11)
 end
 
 ## test large arguments

test/bessely_test.jl

@@ -75,6 +75,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.bessely(BigFloat(v), BigFloat(xx))
     @test isapprox(bessely(v, xx), Float64(sf), rtol=2e-13)
+    @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], Float64(sf), rtol=5e-12)
 end
 
 # test decimal orders
@@ -85,6 +86,7 @@ nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235,
 for v in nu, xx in x
     xx *= v
     @test isapprox(bessely(v, xx), SpecialFunctions.bessely(v, xx), rtol=5e-12)
+    @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], SpecialFunctions.bessely(v, xx), rtol=5e-12)
 end
 
 # need to test accuracy of negative orders and negative arguments and all combinations within