add docs

Author: heltonmc

src/bessely.jl

@@ -181,28 +181,74 @@ function _bessely1_compute(x::Float32)
     end
 end
 
-"""
-    bessely(nu, x::T) where T <: Union{Float32, Float64}
+#              Bessel functions of the second kind of order nu
+#                           bessely(nu, x)
+#
+#    A numerical routine to compute the Bessel function of the second kind Y_{ν}(x) [1]
+#    for real orders and arguments of positive or negative value. The routine is based on several
+#    publications [2, 3, 4, 5] that calculate Y_{ν}(x) for positive arguments and orders where
+#    reflection identities are used to compute negative arguments and orders.
+#
+#    In particular, the reflectance identities for negative integer orders Y_{-n}(x) = (-1)^n * Y_{n}(x) (Eq. 9.1.5; [6])
+#    and for negative noninteger orders Y_{−ν}(x) = cos(πν) * Y_{ν}(x) + sin(πν) * J_{ν}(x) (Eq. 9.1.2; [6]) are used.
+#    For negative arguments of integer order, Y_{n}(-x) = (-1)^n * Y_{n}(x) + (-1)^n * 2im * J_{n}(x) is used and for
+#    noninteger orders, Y_{ν}(−x) = exp(−im*π*ν) * Y_{ν}(x) + 2im * cos(πν) * J_{ν}(x) is used.
+#    For negative orders and arguments the previous identities are combined.
+#
+#    The identities are computed by calling the `bessely_positive_args(nu, x)` function which computes Y_{ν}(x)
+#    for positive arguments and orders. For integer orders up to 250, forward recurrence is used starting from
+#    `bessely0` and `bessely1` routines for calculation of Y_{n}(x) of the zero and first order.
+#    For small arguments, Y_{ν}(x) is calculated from the power series (`bessely_power_series(nu, x`) form of J_{ν}(x) using the connection formula [1].
+#    
+#    When x < ν and ν being reasonably large, the debye asymptotic expansion (Eq. 33; [3]) is used `besseljy_debye(nu, x)`.
+#    When x > ν and x being reasonably large, the Hankel function is calculated from the debye expansion (Eq. 29; [3]) with `hankel_debye(nu, x)`
+#    and Y_{n}(x) is calculated from the imaginary part of the Hankel function.
+#    These expansions are not uniform so are not strictly used when the above inequalities are true, therefore, cutoffs
+#    were determined depending on the desired accuracy. For large arguments x >> ν, the phase functions are used (Eq. 15 [4]) with `besseljy_large_argument(nu, x)`.
+#
+#    For values where the expansions for large arguments and orders are not valid, forward recurrence is employed after shifting the order down
+#    to where these expansions are valid then using recurrence. In general, the routine will be the slowest when ν ≈ x as all methods struggle at this point.
+#    Additionally, the Hankel expansion is only accurate (to Float64 precision) when x > 19 and the power series can only be computed for x < 7 
+#    without loss of precision. Therefore, when x > 19 the Hankel expansion is used to generate starting values after shifting the order down.
+#    When x ∈ (7, 19) and ν ∈ (0, 2) Chebyshev approximation is used with higher orders filled by recurrence.
+#    
+# [1] https://dlmf.nist.gov/10.2#E3
+# [2] Bremer, James. "An algorithm for the rapid numerical evaluation of Bessel functions of real orders and arguments." 
+#     Advances in Computational Mathematics 45.1 (2019): 173-211.
+# [3] Matviyenko, Gregory. "On the evaluation of Bessel functions." 
+#     Applied and Computational Harmonic Analysis 1.1 (1993): 116-135.
+# [4] Heitman, Z., Bremer, J., Rokhlin, V., & Vioreanu, B. (2015). On the asymptotics of Bessel functions in the Fresnel regime. 
+#     Applied and Computational Harmonic Analysis, 39(2), 347-356.
+# [5] Ratis, Yu L., and P. Fernández de Córdoba. "A code to calculate (high order) Bessel functions based on the continued fractions method." 
+#     Computer physics communications 76.3 (1993): 381-388.
+# [6] Abramowitz, Milton, and Irene A. Stegun, eds. Handbook of mathematical functions with formulas, graphs, and mathematical tables. 
+#     Vol. 55. US Government printing office, 1964.
+#
 
-Bessel function of the first kind of order nu, ``J_{nu}(x)``.
-Nu must be real.
+#####
+##### Generic routine for `bessely`
+#####
 """
+    bessely(nu, x::T) where T <: Float64
 
+Bessel function of the first kind of order nu, ``Y_{nu}(x)``.
+nu and x must be real where nu and x can be positive or negative.
+"""
 function bessely(nu::Real, x::T) where T
     isinteger(nu) && return bessely(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
 
     Ynu = bessely_positive_args(abs_nu, abs_x)
-    spi, cpi = sincospi(abs_nu)
     if nu >= zero(T)
         if x >= zero(T)
             return Ynu
         else
-            return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cpi
+            return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cospi(abs_nu)
         end
     else
         Jnu = besselj_positive_args(abs_nu, abs_x)
+        spi, cpi = sincospi(abs_nu)
         if x >= zero(T)
             return Ynu * cpi + Jnu * spi
         else
@@ -224,6 +270,18 @@ function bessely(nu::Integer, x::T) where T
     end
 end
 
+#####
+#####  `bessely` for positive arguments and orders
+#####
+
+"""
+    bessely_positive_args(nu, x::T) where T <: Float64
+
+Bessel function of the first kind of order nu, ``Y_{nu}(x)``.
+nu and x must be real and nu and x must be positive.
+
+No checks on arguments are performed and should only be called if certain nu, x >= 0.
+"""
 function bessely_positive_args(nu, x::T) where T
     iszero(x) && return -T(Inf)
 
@@ -252,16 +310,25 @@ function bessely_positive_args(nu, x::T) where T
     return besselj_up_recurrence(x, imag(hankel_debye(v2, x)), imag(hankel_debye(v2 - 1, x)), v2, nu)[1]
 end
 
+#####
+#####  Power series for Y_{nu}(x)
+#####
+
 # Use power series form of J_v(x) to calculate Y_v(x) with
 # Y_v(x) = (J_v(x)cos(v*π) - J_{-v}(x)) / sin(v*π),    v ~= 0, 1, 2, ...
 # combined to calculate both J_v and J_{-v} in the same loop
 # J_{-v} always converges slower so just check that convergence
 # Combining the loop was faster than using two separate loops
 # 
-# this seems to work well for small arguments x < 7.0 for rel. error ~1e-14
+# this works well for small arguments x < 7.0 for rel. error ~1e-14
 # this also works well for nu > 1.35x - 4.5
 # for nu > 25 more cancellation occurs near integer values
-bessely_series_cutoff(v, x) = (x < 7.0) || v > 1.35*x - 4.5
+"""
+    bessely_power_series(nu, x::T) where T <: Float64
+
+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.
+"""
 function bessely_power_series(v, x::T) where T
     MaxIter = 3000
     out = zero(T)
@@ -281,53 +348,30 @@ function bessely_power_series(v, x::T) where T
     s, c = sincospi(v)
     return (out*c - out2) / s
 end
+bessely_series_cutoff(v, x) = (x < 7.0) || v > 1.35*x - 4.5
 
-#=
-### don't quite have this right
-# probably not needed because we should use debye exapnsion for large nu
-function log_bessely_power_series(v, x::T) where T
-    MaxIter = 2000
-    out = zero(T)
-    out2 = zero(T)
-    a = one(T)
-    b = inv(gamma(1-v))
-    x2 = (x/2)^2
-    for i in 0:MaxIter
-        out += a
-        out2 += b
-        a *= -x2 * inv((i + one(T)) * (v + i + one(T)))
-        b *= -x2 * inv((i + one(T)) * (-v + i + one(T)))
-        (abs(b) < eps(T) * abs(out2)) && break
-    end
-    logout = -loggamma(v + 1) + fma(v, log(x/2), log(out))
-    sign = 1
-
-    if out2 <= zero(T)
-        sign = -1
-        out2 = -out2
-    end
-    logout2 = -v * log(x/2) + log(out2)
-
-    spi, cpi = sincospi(v)
-   
-    #tmp = logout2 + log(abs(sign*(inv(spi)) - exp(logout - logout2) * cpi / spi))
-    tmp = logout + log((-cpi + exp(logout2) - logout) / spi)
+#####
+#####  Chebyshev approximation for Y_{nu}(x)
+#####
 
-    return -exp(tmp)
-end
-=#
+"""
+    bessely_chebyshev(nu, x::T) where T <: Float64
 
-# only implemented for Float64 so far
-besseljy_chebyshev_cutoff(nu, x) = (x <= 19.0 && x >= 6.0)
+Computes ``Y_{nu}(x)`` for medium arguments x ∈ (6, 19) for any positive order using a Chebyshev approximation.
+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]
 end
 
+# only implemented for Float64 so far
+besseljy_chebyshev_cutoff(::Number, x) = (x <= 19.0 && x >= 6.0)
+
 # compute bessely for x ∈ (6, 19) and ν ∈ (0, 2) using chebyshev approximation with a (16, 28) grid (see notes to generate below)
 # optimized to return both (ν, ν + 1) in around the same time, therefore ν must be in (0, 1)
-# no checks are performed on 
+# no checks are performed on arguments
 function bessely_chebyshev_low_orders(v, x)
     # need to rescale inputs according to
     #x0 = (x - lb) * 2 / (ub - lb) - 1
@@ -348,6 +392,7 @@ function clenshaw_chebyshev(x, c)
     end
     return c0 + c1 * x
 end
+
 # to generate the Chebyshev weights
 #=
 using ArbNumerics, FastChebInterp
@@ -376,3 +421,38 @@ const bessely_cheb_weights = (
     (4.9388506677207545e-15, -8.11264980108463e-15, 8.15564695267678e-15, -6.6179014247901285e-15, 1.2560443562777774e-15, -1.1626485837362674e-15, -8.428647433288211e-15, 8.040878625052107e-15, 7.604076442794149e-16, -2.509940602859596e-15, 4.433411320740124e-16, 2.4051574424535735e-16, -7.97263854654364e-17, -8.468775456115063e-18, 3.0977204928150507e-18, -2.763609694292828e-18, -4.528057493615123e-18, 8.090667992714159e-18, -5.7670273029389466e-18, 6.935237082901993e-18, 2.3956539089710175e-18, 3.5297683311975285e-18, -2.6680862870119677e-19, -1.124331067925833e-18, -3.969891564639305e-18, 3.4988194829519685e-18, 8.636411459915624e-19, -1.3895984430563314e-18, 1.7657006809049643e-18),
     (-2.058125495426769e-16, 6.08114692374987e-16, -4.430737266654462e-16, 7.070791826861885e-16, -4.496672312473068e-16, 4.992680927678948e-16, -2.0005398723221536e-16, -4.085152288266179e-16, 2.963046065367326e-16, 1.9081397648228198e-17, -6.065778883208153e-17, 1.2267254817829564e-17, 2.5241246780253446e-18, 5.364001077188067e-19, -1.1151793774136653e-18, -2.647497194530543e-19, -9.99558548618923e-19, 1.9545043306813755e-18, 1.0295729550672243e-18, 5.850363951633064e-19, 2.039938658283451e-19, -1.622964422184252e-18, -2.1240860691038883e-18, 2.0835979355758296e-18, -1.5423943263121515e-19, -5.576729517866199e-19, -3.7259501367076117e-20, -1.8905808347578018e-19, -1.4869058365515489e-18)
 )
+
+#=
+### don't quite have this right (issue with signs)
+# probably not needed because we should use debye exapnsion for large nu
+function log_bessely_power_series(v, x::T) where T
+    MaxIter = 2000
+    out = zero(T)
+    out2 = zero(T)
+    a = one(T)
+    b = inv(gamma(1-v))
+    x2 = (x/2)^2
+    for i in 0:MaxIter
+        out += a
+        out2 += b
+        a *= -x2 * inv((i + one(T)) * (v + i + one(T)))
+        b *= -x2 * inv((i + one(T)) * (-v + i + one(T)))
+        (abs(b) < eps(T) * abs(out2)) && break
+    end
+    logout = -loggamma(v + 1) + fma(v, log(x/2), log(out))
+    sign = 1
+
+    if out2 <= zero(T)
+        sign = -1
+        out2 = -out2
+    end
+    logout2 = -v * log(x/2) + log(out2)
+
+    spi, cpi = sincospi(v)
+   
+    #tmp = logout2 + log(abs(sign*(inv(spi)) - exp(logout - logout2) * cpi / spi))
+    tmp = logout + log((-cpi + exp(logout2) - logout) / spi)
+
+    return -exp(tmp)
+end
+=#