diff --git a/src/besseli.jl b/src/besseli.jl index 714d0e1..c4dc845 100644 --- a/src/besseli.jl +++ b/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 diff --git a/src/besselj.jl b/src/besselj.jl index 74e31ae..afb7b67 100644 --- a/src/besselj.jl +++ b/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 diff --git a/src/besselk.jl b/src/besselk.jl index 7a8caa5..e4bb08d 100644 --- a/src/besselk.jl +++ b/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 diff --git a/src/bessely.jl b/src/bessely.jl index 2412d32..4ee50ca 100644 --- a/src/bessely.jl +++ b/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,7 +245,8 @@ 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) @@ -252,11 +254,11 @@ function bessely(nu::Real, x::T) where T 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) @@ -264,9 +266,19 @@ function bessely(nu::Integer, x::T) where T 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) + end + 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 end diff --git a/test/besseli_test.jl b/test/besseli_test.jl index d25f3dd..72d5f21 100644 --- a/test/besseli_test.jl +++ b/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 diff --git a/test/besselj_test.jl b/test/besselj_test.jl index 238df6e..572bb98 100644 --- a/test/besselj_test.jl +++ b/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) diff --git a/test/besselk_test.jl b/test/besselk_test.jl index a908bad..a1758c9 100644 --- a/test/besselk_test.jl +++ b/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 diff --git a/test/bessely_test.jl b/test/bessely_test.jl index 558c4f6..5a1e9c8 100644 --- a/test/bessely_test.jl +++ b/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)