Added spherical bessel functions (#196)

MasonProtter · web-flow · commit 2820beabb9a4 · 2020-04-05T14:00:21.000-04:00
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "0.10.0"
+version = "0.11.0"
 
 [deps]
 OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
@@ -13,4 +13,4 @@ OpenSpecFun_jll = "0.5.3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test"]
diff --git a/docs/src/functions_list.md b/docs/src/functions_list.md
@@ -32,10 +32,12 @@ SpecialFunctions.besselj0
 SpecialFunctions.besselj1
 SpecialFunctions.besselj
 SpecialFunctions.besseljx
+SpecialFunctions.sphericalbesselj
 SpecialFunctions.bessely0
 SpecialFunctions.bessely1
 SpecialFunctions.bessely
 SpecialFunctions.besselyx
+SpecialFunctions.sphericalbessely
 SpecialFunctions.hankelh1
 SpecialFunctions.hankelh1x
 SpecialFunctions.hankelh2
diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md
@@ -60,10 +60,12 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | [`besselj0(z)`](@ref SpecialFunctions.besselj0)               | `besselj(0,z)`                                                                                                                                                  |
 | [`besselj1(z)`](@ref SpecialFunctions.besselj1)               | `besselj(1,z)`                                                                                                                                                  |
 | [`besseljx(nu,z)`](@ref SpecialFunctions.besseljx)            | scaled Bessel function of the first kind of order `nu` at `z`                                                                                                   |
+| [`sphericalbesselj(nu,z)`](@ref SpecialFunctions.sphericalbesselj) | [Spherical Bessel function](https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn,_yn) of the first kind of order `nu` at `z`       |
 | [`bessely(nu,z)`](@ref SpecialFunctions.bessely)              | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z`                                                        |
 | [`bessely0(z)`](@ref SpecialFunctions.bessely0)               | `bessely(0,z)`                                                                                                                                                  |
 | [`bessely1(z)`](@ref SpecialFunctions.bessely1)               | `bessely(1,z)`                                                                                                                                                  |
 | [`besselyx(nu,z)`](@ref SpecialFunctions.besselyx)            | scaled Bessel function of the second kind of order `nu` at `z`                                                                                                  |
+| [`sphericalbessely(nu,z)`](@ref SpecialFunctions.sphericalbessely) | [Spherical Bessel function](https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn,_yn) of the second kind of order `nu` at `z`      |
 | [`besselh(nu,k,z)`](@ref SpecialFunctions.besselh)            | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the third kind (a.k.a. Hankel function) of order `nu` at `z`; `k` must be either `1` or `2` |
 | [`hankelh1(nu,z)`](@ref SpecialFunctions.hankelh1)            | `besselh(nu, 1, z)`                                                                                                                                             |
 | [`hankelh1x(nu,z)`](@ref SpecialFunctions.hankelh1x)          | scaled `besselh(nu, 1, z)`                                                                                                                                      |
diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl
@@ -16,12 +16,14 @@ export
     besseli,
     besselix,
     besselj,
+    sphericalbesselj,
     besselj0,
     besselj1,
     besseljx,
     besselk,
     besselkx,
     bessely,
+    sphericalbessely,
     bessely0,
     bessely1,
     besselyx,
diff --git a/src/bessel.jl b/src/bessel.jl
@@ -691,6 +691,30 @@ function bessely(n::Integer, x::BigFloat)
     return z
 end
 
+"""
+    sphericalbesselj(nu, x)
+
+Spherical bessel function of the first kind at order `nu`, ``j_ν(x)``. This is the non-singular
+solution to the radial part of the Helmholz equation in spherical coordinates.
+"""
+function sphericalbesselj(nu, x::T) where {T}
+    besselj_nuhalf_x = besselj(nu + one(nu)/2, x)
+    if abs(x) ≤ sqrt(eps(real(zero(besselj_nuhalf_x))))
+        nu == 0 ? one(besselj_nuhalf_x) : zero(besselj_nuhalf_x)
+    else
+        √((float(T))(π)/2x) * besselj_nuhalf_x
+    end
+end
+
+"""
+    sphericalbessely(nu, x)
+
+Spherical bessel function of the second kind at order `nu`, ``y_ν(x)``. This is the singular
+solution to the radial part of the Helmholz equation in spherical coordinates. Sometimes
+known as a spherical Neumann function.
+"""
+sphericalbessely(nu, x::T) where {T} = √((float(T))(π)/2x) * bessely(nu + 1//2, x)
+
 """
     hankelh1(nu, x)
 
diff --git a/test/bessel.jl b/test/bessel.jl
@@ -189,6 +189,33 @@ end
         end
     end
 
+    @testset "sphericalbesselj" begin
+        @test sphericalbesselj(1, 1)      ≈ 0.3011686789397568
+        @test sphericalbesselj(10, 5.5)   ≈ 0.0009369210263385842
+        @test sphericalbesselj(1.25, 5.5) ≈ -0.1123558799930763
+        @test sphericalbesselj(1.25, -5.5+0im) ≈ 0.079447604649286 + 0.079447604649286im
+
+        @test sphericalbesselj(0, 0.01) ≈ 0.999983333416666
+        @test sphericalbesselj(0, 0)    == 1.0
+        @test sphericalbesselj(1, 0)    == 0.0
+        @test sphericalbesselj(1, 0.01) ≈ 0.003333300000119047
+        
+        @test_throws DomainError sphericalbesselj(1.25, -5.5)
+    end
+
+    @testset "sphericalbessely" begin
+        @test sphericalbessely(1, 1)      ≈ -1.381773290676036
+        @test sphericalbessely(10, 5.5)   ≈ -10.89087037026398
+        @test sphericalbessely(1.25, 5.5) ≈ 0.148322390312228
+        @test sphericalbessely(1.25, -5.5+0im) ≈ -0.054015441306998 - 0.104879767991574im
+
+        @test sphericalbessely(0, 1e-5) ≈ -99999.9999950000000
+        @test sphericalbessely(1, 1e-5) ≈ -1e10
+        
+        @test_throws DomainError sphericalbessely(1.25, -5.5)
+        @test_throws AmosException sphericalbessely(1, 0)
+    end
+
     @testset "besselhx" begin
         for elty in [Complex{Float16},Complex{Float32},Complex{Float64}]
             z = convert(elty, 1.0 + 1.9im)
