add hankel functions

Author: heltonmc

src/hankel.jl

@@ -194,6 +194,25 @@ function besselh(nu::Real, k::Integer, x)
     end
 end
 
+function besselh(nu::AbstractRange, k::Integer, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    if nu[end] < x
+        out = Vector{Complex{T}}(undef, length(nu))
+        out[1], out[2] = besselh(nu[1], k, x), besselh(nu[2], k, x)
+        return besselj_up_recurrence!(out, x, nu)
+    else
+        Jn = besselj(nu, x)
+        Yn = bessely(nu, x)
+        if k == 1
+            return complex.(Jn, Yn)
+        elseif k == 2
+            return complex.(Jn, -Yn)
+        else
+            throw(ArgumentError("k must be 1 or 2"))
+        end
+    end
+end
+
 """
     hankelh1(nu, x)
 Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``.

test/hankel_test.jl

@@ -26,3 +26,8 @@ v, x = 14.3, 29.4
 @test isapprox(hankelh2(v, x), SpecialFunctions.hankelh2(v, x), rtol=2e-13)
 @test isapprox(besselh(v, 1, x), SpecialFunctions.besselh(v, 1, x), rtol=2e-13)
 @test isapprox(besselh(v, 2, x), SpecialFunctions.besselh(v, 2, x), rtol=2e-13)
+
+@test isapprox(hankelh1(1:50, 10.0), SpecialFunctions.hankelh1.(1:50, 10.0), rtol=2e-13)
+@test isapprox(hankelh1(0.5:1:25.5, 15.0), SpecialFunctions.hankelh1.(0.5:1:25.5, 15.0), rtol=2e-13)
+@test isapprox(hankelh1(1:50, 100.0), SpecialFunctions.hankelh1.(1:50, 100.0), rtol=2e-13)
+@test isapprox(hankelh2(1:50, 10.0), SpecialFunctions.hankelh2.(1:50, 10.0), rtol=2e-13)