From 30e0b3dd7c5adcfe9478ee6b042263499ddbcee1 Mon Sep 17 00:00:00 2001
From: jondea <jondea@thedeakin.net>
Date: Mon, 12 Apr 2021 00:23:27 +0100
Subject: [PATCH] Implement invhankelh1 defined by passing through and tests
 for it

---
 src/InverseHankelFunction.jl |  8 ++++++-
 src/passingthrough.jl        | 32 +++++++++++----------------
 test/runtests.jl             | 42 +++++++++++++++++++++++++++++++++++-
 3 files changed, 61 insertions(+), 21 deletions(-)

diff --git a/src/InverseHankelFunction.jl b/src/InverseHankelFunction.jl
index 2bb5ebc..3896a0f 100644
--- a/src/InverseHankelFunction.jl
+++ b/src/InverseHankelFunction.jl
@@ -3,7 +3,11 @@ module InverseHankelFunction
 import SpecialFunctions: hankelh1, gamma
 import LambertW: lambertw
 
-export invhankelh1, PassingThrough
+export invhankelh1,
+       PassingThrough,
+       isfinite,
+       iszero,
+       isnan
 
 export hankelh1n,
        HankelH1N
@@ -34,4 +38,6 @@ include("hankel_asymptotics.jl")
 
 include("invhankel.jl")
 
+include("passingthrough.jl")
+
 end # module
diff --git a/src/passingthrough.jl b/src/passingthrough.jl
index dbe6078..01bc847 100644
--- a/src/passingthrough.jl
+++ b/src/passingthrough.jl
@@ -3,9 +3,14 @@ struct PassingThrough{T}
     point::T
 end
 
+import Base.isfinite
 isfinite(p::PassingThrough) = isfinite(p.point)
+
+import Base.iszero
 iszero(p::PassingThrough) = iszero(p.point)
 
+import Base.isnan
+isnan(p::PassingThrough) = isnan(p.point)
 
 """
     invhankelh1(ν, h, p::PassingThrough)
@@ -15,29 +20,18 @@ Find z such that H^{(1)}_\\nu(z) = h and the solution branch passes through p, t
 function invhankelh1(ν::Integer, h::Number, p::PassingThrough)
 
     # Handle some special cases
-    if isinf(p)
-        throw(DomainError(p.point, "Multiple branches of invhankelh1 pass through $(p.point), find another value that this branch passes through"))
+    if isnan(p)
+        throw(DomainError(p.point, "invhankelh1 never passes through $(p.point) because it is not a number"))
     end
-    if iszero(p)
-        throw(DomainError(p.point, "Multiple branches of invhankelh1 pass through $(p.point), find another value that this branch passes through"))
+    if !isfinite(p)
+        throw(DomainError(p.point, "Multiple branches of invhankelh1 pass through non finite ($(p.point)) find another value that this branch passes through"))
     end
-    if isnan(p)
-        throw(DomainError(p.point, "invhankelh1 never passes through $(p.point)"))
+    if iszero(p)
+        throw(DomainError(p.point, "Multiple branches of invhankelh1 pass through zero ($(p.point)) find another value that this branch passes through"))
     end
 
     # Starting point for our continuation
-    z = p.point
-
-    SMALL_ARGUMENT_THRESHOLD = 0.5
-    LARGE_ARGUMENT_THRESHOLD = 2
-    if abs(z)/sqrt(abs(ν)+1) < SMALL_ARGUMENT_THRESHOLD
-
-    elseif abs(z)/sqrt(abs(ν)+1) > LARGE_ARGUMENT_THRESHOLD
-
-    end
-
-
-
-    return 1.0
+    z₀ = p.point
 
+    return invhankelh1n(ν, z₀, h/hankelh1(ν,z₀))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 62f9c52..b5bc224 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -107,7 +107,7 @@ include("utils.jl")
 
 end
 
-@testset "Normalised Hankel function" begin
+@testset "Normalised Hankel function        " begin
     ν = 1
     z₀ = 2.0 - 0.3im
     z = 2.0 + 0.1im
@@ -239,6 +239,46 @@ end
 
 end
 
+@testset "Passing through                   " begin
+    for ν in [0, 1, 4]
+        for p in PassingThrough.([1.0, 2.0, 4.0, 1.0+0.5im, 2.0-0.2im])
+            h(z) = hankelh1(ν, z)
+            h⁻¹(z) = invhankelh1(ν, z, p)
+
+            # Check that invhankelh1 does indeed pass through this point
+            @test h⁻¹(h(p.point)) ≈ p.point
+
+            for hz in [0.9, 0.34, 0.01, 1.4, 1.0-0.5im]
+                h⁻¹hz = h⁻¹(hz)
+                @test h(h⁻¹hz) ≈ hz
+            end
+        end
+    end
+
+    @testset "Domain errors" begin
+        ν = 4
+        h = 3.4 + 0.1im
+
+        # NaNs
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN - 1.0*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN - Inf*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN - NaN*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(10  - NaN*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN))
+
+        # Any kind of infinity
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(Inf + 0im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(1.0 + Inf*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(Inf + -Inf*im))
+
+        # Zeros
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough( 0.0 + 0.0*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough( 0.0 - 0.0*im))
+        @test_throws DomainError invhankelh1(ν, h, PassingThrough(-0.0 - 0.0*im))
+    end
+
+end
+
 @testset "Derivatives of Hankel functions   " begin
 
     for ν in [0, 1, 4, 10]