Skip to content

Commit a3c767f

Browse files
heltonmcheltonmc
committed
add sphericalbesseli
1 parent 43b5e20 commit a3c767f

File tree

3 files changed

+94
-8
lines changed

3 files changed

+94
-8
lines changed

‎src/besseli.jl

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -170,6 +170,10 @@ _besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x)))
170170

171171
function _besseli(nu, x::T) where T <: Union{Float32, Float64}
172172
isinteger(nu) && return _besseli(Int(nu), x)
173+
if ~isfinite(x)
174+
isnan(x) && return x
175+
isinf(x) && return x
176+
end
173177
abs_nu = abs(nu)
174178
abs_x = abs(x)
175179

@@ -192,6 +196,10 @@ function _besseli(nu, x::T) where T <: Union{Float32, Float64}
192196
end
193197
end
194198
function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64}
199+
if ~isfinite(x)
200+
isnan(x) && return x
201+
isinf(x) && return x
202+
end
195203
abs_nu = abs(nu)
196204
abs_x = abs(x)
197205
sg = iseven(abs_nu) ? 1 : -1
@@ -211,7 +219,6 @@ Modified Bessel function of the first kind of order nu, ``I_{nu}(x)`` for positi
211219
function besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64}
212220
iszero(nu) && return besseli0(x)
213221
isone(nu) && return besseli1(x)
214-
isinf(x) && return T(Inf)
215222

216223
# use large argument expansion if x >> nu
217224
besseli_large_argument_cutoff(nu, x) && return besseli_large_argument(nu, x)

‎src/modifiedsphericalbessel.jl

Lines changed: 48 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,16 @@
11
# Modified Spherical Bessel functions
22
#
3-
# sphericalbesselk(nu, x)
3+
# sphericalbesseli(nu, x), sphericalbesselk(nu, x)
44
#
5-
# A numerical routine to compute the modified spherical bessel functions of the second kind.
6-
# For moderate sized integer orders, forward recurrence is used starting from explicit formulas for k0(x) = exp(-x) / x and k1(x) = k0(x) * (x+1) / x [1].
5+
# A numerical routine to compute the modified spherical bessel functions of the first and second kind.
6+
# The modified spherical bessel function of the first kind is computed using the power series for small arguments,
7+
# explicit formulas for (nu=0,1,2), and using its relation to besseli for other arguments [1].
8+
# The modified bessel function of the second kind is computed for small to moderate integer orders using forward recurrence starting from explicit formulas for k0(x) = exp(-x) / x and k1(x) = k0(x) * (x+1) / x [2].
79
# Large orders are determined from the uniform asymptotic expansions (see src/besselk.jl for details)
8-
# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [1].
9-
#
10-
# [1] https://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html
10+
# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [2].
11+
#
12+
# [1] https://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html
13+
# [2] https://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html
1114
#
1215
"""
1316
sphericalbesselk(nu, x::T) where T <: {Float32, Float64}
@@ -54,3 +57,42 @@ function sphericalbesselk_int(v::Int, x)
5457
end
5558
b1
5659
end
60+
61+
"""
62+
sphericalbesseli(nu, x::T) where T <: {Float32, Float64}
63+
64+
Computes `i_{ν}(x)`, the modified first-kind spherical Bessel function.
65+
"""
66+
sphericalbesseli(nu::Real, x::Real) = _sphericalbesseli(nu, float(x))
67+
68+
_sphericalbesseli(nu, x::Float16) = Float16(_sphericalbesseli(nu, Float32(x)))
69+
70+
function _sphericalbesseli(nu, x::T) where T <: Union{Float32, Float64}
71+
isinf(x) && return x
72+
x < zero(x) && throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
73+
74+
sphericalbesselj_small_args_cutoff(nu, x::T) && return sphericalbesseli_small_args(nu, x)
75+
isinteger(nu) && return _sphericalbesseli_small_orders(Int(nu), x)
76+
return SQPIO2(T)*besseli(nu+1/2, x) / sqrt(x)
77+
end
78+
79+
function _sphericalbesseli_small_orders(nu::Integer, x::T) where T
80+
# prone to cancellation in the subtraction
81+
# best to expand and group
82+
nu_abs = abs(nu)
83+
x2 = x*x
84+
sinhx = sinh(x)
85+
coshx = cosh(x)
86+
nu_abs == 0 && return sinhx / x
87+
nu_abs == 1 && return (x*coshx - sinhx) / x2
88+
nu_abs == 2 && return (x2*sinhx + 3*(sinhx - x*coshx)) / (x2*x)
89+
return SQPIO2(T)*besseli(nu+1/2, x) / sqrt(x)
90+
end
91+
92+
function sphericalbesseli_small_args(nu, x::T) where T
93+
iszero(x) && return iszero(nu) ? one(T) : x
94+
x2 = x^2 / 4
95+
coef = evalpoly(x2, (1, inv(T(3)/2 + nu), inv(5 + nu), inv(T(21)/2 + nu), inv(18 + nu)))
96+
a = SQPIO2(T) / (gamma(T(3)/2 + nu) * 2^(nu + T(1)/2))
97+
return x^nu * a * coef
98+
end

‎test/sphericalbessel_test.jl

Lines changed: 38 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,14 @@ x = 1e-15
1111
@test Bessels.sphericalbesselk(5.5, x) ≈ SpecialFunctions.besselk(5.5 + 1/2, x) * sqrt( 2 / (x*pi))
1212
@test Bessels.sphericalbesselk(10, x) ≈ SpecialFunctions.besselk(10 + 1/2, x) * sqrt( 2 / (x*pi))
1313

14+
x = 1e-20
15+
@test Bessels.sphericalbesseli(0, x) ≈ SpecialFunctions.besseli(0 + 1/2, x) * sqrt( pi / (x*2))
16+
@test Bessels.sphericalbesseli(1, x) ≈ SpecialFunctions.besseli(1 + 1/2, x) * sqrt( pi / (x*2))
17+
@test Bessels.sphericalbesseli(2, x) ≈ SpecialFunctions.besseli(2 + 1/2, x) * sqrt( pi / (x*2))
18+
@test Bessels.sphericalbesseli(3, x) ≈ SpecialFunctions.besseli(3 + 1/2, x) * sqrt( pi / (x*2))
19+
@test Bessels.sphericalbesseli(4, x) ≈ SpecialFunctions.besseli(4 + 1/2, x) * sqrt( pi / (x*2))
20+
@test Bessels.sphericalbesseli(6.5, x) ≈ SpecialFunctions.besseli(6.5 + 1/2, x) * sqrt( pi / (x*2))
21+
1422
# test zero
1523
@test isone(Bessels.sphericalbesselj(0, 0.0))
1624
@test iszero(Bessels.sphericalbesselj(3, 0.0))
@@ -26,6 +34,14 @@ x = 1e-15
2634
@test isinf(Bessels.sphericalbesselk(4, 0.0))
2735
@test isinf(Bessels.sphericalbesselk(10.2, 0.0))
2836

37+
x = 0.0
38+
@test isone(Bessels.sphericalbesseli(0, x))
39+
@test iszero(Bessels.sphericalbesseli(1, x))
40+
@test iszero(Bessels.sphericalbesseli(2, x))
41+
@test iszero(Bessels.sphericalbesseli(3, x))
42+
@test iszero(Bessels.sphericalbesseli(4, x))
43+
@test iszero(Bessels.sphericalbesseli(6.4, x))
44+
2945
# test Inf
3046
@test iszero(Bessels.sphericalbesselj(1, Inf))
3147
@test iszero(Bessels.sphericalbesselj(10.2, Inf))
@@ -36,6 +52,14 @@ x = 1e-15
3652
@test iszero(Bessels.sphericalbesselk(4, Inf))
3753
@test iszero(Bessels.sphericalbesselk(10.2, Inf))
3854

55+
x = Inf
56+
@test isinf(Bessels.sphericalbesseli(0, x))
57+
@test isinf(Bessels.sphericalbesseli(1, x))
58+
@test isinf(Bessels.sphericalbesseli(2, x))
59+
@test isinf(Bessels.sphericalbesseli(3, x))
60+
@test isinf(Bessels.sphericalbesseli(4, x))
61+
@test isinf(Bessels.sphericalbesseli(6.4, x))
62+
3963
# test NaN
4064
@test isnan(Bessels.sphericalbesselj(1.4, NaN))
4165
@test isnan(Bessels.sphericalbesselj(4.0, NaN))
@@ -45,6 +69,14 @@ x = 1e-15
4569
@test isnan(Bessels.sphericalbesselk(1.4, NaN))
4670
@test isnan(Bessels.sphericalbesselk(4.0, NaN))
4771

72+
x = NaN
73+
@test isnan(Bessels.sphericalbesseli(0, x))
74+
@test isnan(Bessels.sphericalbesseli(1, x))
75+
@test isnan(Bessels.sphericalbesseli(2, x))
76+
@test isnan(Bessels.sphericalbesseli(3, x))
77+
@test isnan(Bessels.sphericalbesseli(4, x))
78+
@test isnan(Bessels.sphericalbesseli(6.4, x))
79+
4880
# test Float16, Float32 types
4981
@test Bessels.sphericalbesselj(Float16(1.4), Float16(1.2)) isa Float16
5082
@test Bessels.sphericalbessely(Float16(1.4), Float16(1.2)) isa Float16
@@ -54,16 +86,21 @@ x = 1e-15
5486
@test Bessels.sphericalbesselk(Float16(1.4), Float16(1.2)) isa Float16
5587
@test Bessels.sphericalbesselk(1.0f0, 1.2f0) isa Float32
5688

57-
for x in 0.5:1.0:100.0, v in [0, 1, 5.5, 8.2, 10]
89+
@test Bessels.sphericalbesseli(Float16(1.4), Float16(1.2)) isa Float16
90+
@test Bessels.sphericalbesseli(1.0f0, 1.2f0) isa Float32
91+
92+
for x in 0.5:1.5:100.0, v in [0, 1, 2, 3, 4, 5.5, 8.2, 10]
5893
@test isapprox(Bessels.sphericalbesselj(v, x), SpecialFunctions.sphericalbesselj(v, x), rtol=1e-12)
5994
@test isapprox(Bessels.sphericalbessely(v, x), SpecialFunctions.sphericalbessely(v, x), rtol=1e-12)
6095
@test isapprox(Bessels.sphericalbesselk(v, x), SpecialFunctions.besselk(v+1/2, x) * sqrt( 2 / (x*pi)), rtol=1e-12)
96+
@test isapprox(Bessels.sphericalbesseli(v, x), SpecialFunctions.besseli(v+1/2, x) * sqrt( pi / (x*2)), rtol=1e-12)
6197
end
6298

6399
for x in 5.5:4.0:160.0, v in [20, 25.0, 32.4, 40.0, 45.12, 50.0, 55.2, 60.124, 70.23, 75.0, 80.0, 92.3, 100.0, 120.0]
64100
@test isapprox(Bessels.sphericalbesselj(v, x), SpecialFunctions.sphericalbesselj(v, x), rtol=3e-12)
65101
@test isapprox(Bessels.sphericalbessely(v, x), SpecialFunctions.sphericalbessely(v, x), rtol=3e-12)
66102
@test isapprox(Bessels.sphericalbesselk(v, x), SpecialFunctions.besselk(v+1/2, x) * sqrt( 2 / (x*pi)), rtol=1e-12)
103+
@test isapprox(Bessels.sphericalbesseli(v, x), SpecialFunctions.besseli(v+1/2, x) * sqrt( pi / (x*2)), rtol=1e-12)
67104
end
68105

69106
@test isapprox(Bessels.sphericalbessely(270, 240.0), SpecialFunctions.sphericalbessely(270, 240.0), rtol=3e-12)

0 commit comments

Comments
 (0)