remove complex output for real inputs

Author: heltonmc

src/besseli.jl

@@ -173,15 +173,17 @@ function besseli(nu::Real, x::T) where T
         if x >= 0
             return besseli_positive_args(abs_nu, abs_x)
         else
-            return cispi(abs_nu) * besseli_positive_args(abs_nu, abs_x)
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return cispi(abs_nu) * besseli_positive_args(abs_nu, abs_x)
         end
     else
         if x >= 0
             return besseli_positive_args(abs_nu, abs_x) + 2 / π * sinpi(abs_nu) * besselk_positive_args(abs_nu, abs_x)
         else
-            Iv = besseli_positive_args(abs_nu, abs_x)
-            Kv = besselk_positive_args(abs_nu, abs_x)
-            return cispi(abs_nu) * Iv + 2 / π * sinpi(abs_nu) * (cispi(-abs_nu) * Kv - im * π * Iv)
+            #Iv = besseli_positive_args(abs_nu, abs_x)
+            #Kv = besselk_positive_args(abs_nu, abs_x)
+            #return cispi(abs_nu) * Iv + 2 / π * sinpi(abs_nu) * (cispi(-abs_nu) * Kv - im * π * Iv)
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
         end
     end
 end

src/besselj.jl

@@ -199,12 +199,22 @@ function besselj(nu::Real, x::T) where T
 
     Jnu = besselj_positive_args(abs_nu, abs_x)
     if nu >= zero(T)
-        return x >= zero(T) ? Jnu : Jnu * cispi(abs_nu)
+        if x >= zero(T)
+            return Jnu
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Jnu * cispi(abs_nu)
+        end
     else
         Ynu = bessely_positive_args(abs_nu, abs_x)
         spi, cpi = sincospi(abs_nu)
         out = Jnu * cpi - Ynu * spi
-        return x >= zero(T) ? out : out * cispi(nu)
+        if x >= zero(T)
+            return out
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return out * cispi(nu)
+        end
     end
 end
 

src/besselk.jl

@@ -193,7 +193,8 @@ function besselk(nu::Real, x::T) where T
     if x >= 0
         return besselk_positive_args(abs_nu, abs_x)
     else
-        return cispi(-abs_nu)*besselk_positive_args(abs_nu, abs_x) - besseli_positive_args(abs_nu, abs_x) * im * π
+        return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+        #return cispi(-abs_nu)*besselk_positive_args(abs_nu, abs_x) - besseli_positive_args(abs_nu, abs_x) * im * π
     end
 end
 function besselk(nu::Integer, x::T) where T
@@ -204,7 +205,8 @@ function besselk(nu::Integer, x::T) where T
     if x >= 0
         return besselk_positive_args(abs_nu, abs_x)
     else
-        return sg * besselk_positive_args(abs_nu, abs_x) - im * π * besseli_positive_args(abs_nu, abs_x)
+        return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+        #return sg * besselk_positive_args(abs_nu, abs_x) - im * π * besseli_positive_args(abs_nu, abs_x)
     end
 end
 

src/bessely.jl

@@ -234,6 +234,7 @@ end
 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)
@@ -244,29 +245,38 @@ function bessely(nu::Real, x::T) where T
         if x >= zero(T)
             return Ynu
         else
-            return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cospi(abs_nu)
+            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
         Jnu = besselj_positive_args(abs_nu, abs_x)
         spi, cpi = sincospi(abs_nu)
         if x >= zero(T)
             return Ynu * cpi + Jnu * spi
         else
-            return cpi * (Ynu * cispi(nu) + 2im * Jnu * cpi) + Jnu * spi * cispi(abs_nu)
+            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 bessely(nu::Integer, x::T) where T
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
 
     Ynu = bessely_positive_args(abs_nu, abs_x)
     if nu >= zero(T)
-        return x >= zero(T) ? Ynu : Ynu * sg + 2im * sg * besselj_positive_args(abs_nu, abs_x)
-    elseif nu < zero(T)
-        return x >= zero(T) ? Ynu * sg : Ynu + 2im * besselj_positive_args(abs_nu, abs_x)
+        if x >= zero(T)
+            return Ynu
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Ynu * sg + 2im * sg * besselj_positive_args(abs_nu, abs_x)
+    else
+        if x >= zero(T)
+            return Ynu * sg
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Ynu + 2im * besselj_positive_args(abs_nu, abs_x)
     end
 end
 

test/besseli_test.jl

@@ -111,7 +111,7 @@ end
 @test besseli(v,x) ≈ -1.995631678207200756444e-14
 
 (v,x) = 12.6, -3.0
-@test besseli(v,x) ≈ -2.725684975265100582482e-8 + 8.388795775899337839603e-8 * im
+#@test besseli(v,x) ≈ -2.725684975265100582482e-8 + 8.388795775899337839603e-8 * im
 
 (v, x) = -8.0, 4.2
 @test besseli(v,x) ≈ 0.0151395115677545706602449919627
@@ -126,4 +126,4 @@ end
 @test besseli(v,x) ≈ 0.2892290867115615816280234648
 
 (v, x) = -14.6, -10.6
-@test besseli(v,x) ≈ -0.157582642056898598750175404443 - 0.484989503203097528858271270828*im
+#@test besseli(v,x) ≈ -0.157582642056898598750175404443 - 0.484989503203097528858271270828*im

test/besselj_test.jl

@@ -118,7 +118,7 @@ end
 # values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica
 # need to also account for different branches when nu isa integer
 nu = -9.102; x = -12.48
-@test isapprox(besselj(nu, x), 0.09842356047575545808128 -0.03266486217437818486161im, rtol=8e-14)
+#@test isapprox(besselj(nu, x), 0.09842356047575545808128 -0.03266486217437818486161im, rtol=8e-14)
 nu = -5.0; x = -5.1
 @test isapprox(besselj(nu, x), 0.2740038554704588327387, rtol=8e-14)
 nu = -7.3; x = 19.1
@@ -128,4 +128,4 @@ nu = -14.0; x = 21.3
 nu = 13.0; x = -8.5
 @test isapprox(besselj(nu, x), -0.006128034621313167000171, rtol=8e-14)
 nu = 17.45; x = -16.23
-@test isapprox(besselj(nu, x), -0.01607335977752705869797 -0.1014831996412783806255im, rtol=8e-14)
+#@test isapprox(besselj(nu, x), -0.01607335977752705869797 -0.1014831996412783806255im, rtol=8e-14)

test/besselk_test.jl

@@ -126,10 +126,10 @@ end
 @test besselk(v,x) ≈ 56331.504348755621996013084096
 
 (v,x) = 13.0, -1.0
-@test besselk(v,x) ≈ -1921576392792.994084565 - 6.26946181952681217841e-14*im
+#@test besselk(v,x) ≈ -1921576392792.994084565 - 6.26946181952681217841e-14*im
 
 (v,x) = 12.6, -3.0
-@test besselk(v,x) ≈ -135222.7926354826692727 - 416172.9627453652120223*im
+#@test besselk(v,x) ≈ -135222.7926354826692727 - 416172.9627453652120223*im
 
 (v, x) = -8.0, 4.2
 @test besselk(v,x) ≈ 3.65165949039881135495282699061
@@ -141,7 +141,7 @@ end
 @test besselk(v,x) ≈ 0.299085139926840649406079315812
 
 (v, x) = -14.0, -9.9
-@test besselk(v,x) ≈ 0.100786833375605803570325345603 - 0.90863997401752715470886289641*im
+#@test besselk(v,x) ≈ 0.100786833375605803570325345603 - 0.90863997401752715470886289641*im
 
 (v, x) = -14.6, -10.6
-@test besselk(v,x) ≈ -0.0180385087581148387140033906859 - 1.54653251445680014758965158559*im
+#@test besselk(v,x) ≈ -0.0180385087581148387140033906859 - 1.54653251445680014758965158559*im

test/bessely_test.jl

@@ -92,14 +92,14 @@ end
 # values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica
 # need to also account for different branches when nu isa integer
 nu = -2.3; x = -2.4
-@test isapprox(bessely(nu, x), -0.0179769671833457636186 + 0.7394120337538928700168im, rtol=5e-14)
+#@test isapprox(bessely(nu, x), -0.0179769671833457636186 + 0.7394120337538928700168im, rtol=5e-14)
 nu = -4.0; x = -12.6
-@test isapprox(bessely(nu, x), -0.02845106816742465563357 + 0.4577922229605476882792im, rtol=5e-14)
+#@test isapprox(bessely(nu, x), -0.02845106816742465563357 + 0.4577922229605476882792im, rtol=5e-14)
 nu = -6.2; x = 18.6
 @test isapprox(bessely(nu, x), -0.05880321550673669650027, rtol=5e-14)
 nu = -8.0; x = 23.2
 @test isapprox(bessely(nu, x),-0.166071277370329242677, rtol=5e-14)
 nu = 11.0; x = -8.2
-@test isapprox(bessely(nu, x),1.438494049708244558901 -0.06222860017222637350092im, rtol=5e-14)
+#@test isapprox(bessely(nu, x),1.438494049708244558901 -0.06222860017222637350092im, rtol=5e-14)
 nu = 13.678; x = -12.98
-@test isapprox(bessely(nu, x),-0.2227392320508850571009 -0.2085585256158188848322im, rtol=5e-14)
+#@test isapprox(bessely(nu, x),-0.2227392320508850571009 -0.2085585256158188848322im, rtol=5e-14)