add more docs and udpate compat

Project.toml

@@ -3,11 +3,8 @@ uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 authors = ["Michael Helton <[email protected]> and contributors"]
 version = "0.1.0"
 
-[deps]
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-
 [compat]
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/Bessels.jl

@@ -1,7 +1,5 @@
 module Bessels
 
-using SpecialFunctions: loggamma
-
 export besselj0
 export besselj1
 export besselj

src/U_polynomials.jl

@@ -1,5 +1,19 @@
-besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && x > 15
-# valid when x < v (uniform asymptotic expansions)
+#                 Debye asymptotic expansions
+#                `besseljy_debye`, `hankel_debye`
+#
+#  This file contains the debye asymptotic asymptotic expansions for large orders.
+#  These routines can be used to calculate `besselj`, `bessely`, `besselk`, `besseli`
+#  and the Hankel functions. These are uniform expansions for `besselk` and `besseli` for large orders.
+#  The forms used for `besselj` and `bessely` as well as the Hankel functions follow the notation by Matviyenko [1]
+#  but also see similar forms provided by NIST 10.19. All of these routines use a core routine for calculating
+#  the U-polynomials (NIST 10.41.E9) where the coefficients are dependent on each function. The forms for the modified
+#  Bessel functions `besselk` and `besseli` follow NIST 10.41 [3].
+#
+# [1] Matviyenko, Gregory. "On the evaluation of Bessel functions." 
+#     Applied and Computational Harmonic Analysis 1.1 (1993): 116-135.
+# [2] http://dlmf.nist.gov/10.41.E9
+# [3] https://dlmf.nist.gov/10.41
+
 """
     besseljy_debye(nu, x::T)
 
@@ -26,9 +40,9 @@ function besseljy_debye(v, x)
 
     return coef_Jn * Uk_Jn, coef_Yn * Uk_Yn
 end
-hankel_debye_cutoff(nu, x) = nu < 0.2 + x + Base.Math._approx_cbrt(-411*x)
 
-# valid when v < x (uniform asymptotic expansions)
+besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && x > 15
+
 """
     hankel_debye(nu, x::T)
 
@@ -54,6 +68,8 @@ function hankel_debye(v, x::T) where T
     return coef_Yn * Uk_Yn
 end
 
+hankel_debye_cutoff(nu, x) = nu < 0.2 + x + Base.Math._approx_cbrt(-411*x)
+
 function Uk_poly_Jn(p, v, p2, x::T) where T <: Float64
     if v > 5.0 + 1.00033*x + Base.Math._approx_cbrt(1427.61*x)
         return Uk_poly10(p, v, p2)

src/asymptotics.jl

@@ -1,3 +1,13 @@
+#             Asymptotics expansions for x > nu
+#                `besseljy_large_argument`
+#
+#  This file contains the asymptotic asymptotic expansions for large arguments 
+#  which is accurate in the fresnel regime |x| > nu following the derivations given by Heitman [1].
+#  These routines can be used to calculate `besselj` and `bessely` using phase functions.
+#
+# [1] 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.
+
 #besseljy_large_argument_min(::Type{Float32}) = 15.0f0
 besseljy_large_argument_min(::Type{Float64}) = 20.0
 besseljy_large_argument_min(::Type{T}) where T <: AbstractFloat = 40.0
@@ -6,7 +16,12 @@ besseljy_large_argument_min(::Type{T}) where T <: AbstractFloat = 40.0
 besseljy_large_argument_cutoff(v, x::Float64) = (x > 1.65*v && x > besseljy_large_argument_min(Float64))
 besseljy_large_argument_cutoff(v, x::T) where T = (x > 4*v && x > besseljy_large_argument_min(T))
 
+"""
+    besseljy_large_argument(nu, x::T)
 
+Asymptotic expansions for large arguments valid when x > 1.6*nu and x > 20.0.
+Returns both (besselj(nu, x), bessely(nu, x)).
+"""
 function besseljy_large_argument(v, x::T) where T
     # gives both (besselj, bessely) for x > 1.6*v
     α, αp = _α_αp_asymptotic(v, x)