1# 2# Tests of spherical Bessel functions. 3# 4import numpy as np 5from numpy.testing import (assert_almost_equal, assert_allclose, 6 assert_array_almost_equal, suppress_warnings) 7import pytest 8from numpy import sin, cos, sinh, cosh, exp, inf, nan, r_, pi 9 10from scipy.special import spherical_jn, spherical_yn, spherical_in, spherical_kn 11from scipy.integrate import quad 12 13 14class TestSphericalJn: 15 def test_spherical_jn_exact(self): 16 # https://dlmf.nist.gov/10.49.E3 17 # Note: exact expression is numerically stable only for small 18 # n or z >> n. 19 x = np.array([0.12, 1.23, 12.34, 123.45, 1234.5]) 20 assert_allclose(spherical_jn(2, x), 21 (-1/x + 3/x**3)*sin(x) - 3/x**2*cos(x)) 22 23 def test_spherical_jn_recurrence_complex(self): 24 # https://dlmf.nist.gov/10.51.E1 25 n = np.array([1, 2, 3, 7, 12]) 26 x = 1.1 + 1.5j 27 assert_allclose(spherical_jn(n - 1, x) + spherical_jn(n + 1, x), 28 (2*n + 1)/x*spherical_jn(n, x)) 29 30 def test_spherical_jn_recurrence_real(self): 31 # https://dlmf.nist.gov/10.51.E1 32 n = np.array([1, 2, 3, 7, 12]) 33 x = 0.12 34 assert_allclose(spherical_jn(n - 1, x) + spherical_jn(n + 1,x), 35 (2*n + 1)/x*spherical_jn(n, x)) 36 37 def test_spherical_jn_inf_real(self): 38 # https://dlmf.nist.gov/10.52.E3 39 n = 6 40 x = np.array([-inf, inf]) 41 assert_allclose(spherical_jn(n, x), np.array([0, 0])) 42 43 def test_spherical_jn_inf_complex(self): 44 # https://dlmf.nist.gov/10.52.E3 45 n = 7 46 x = np.array([-inf + 0j, inf + 0j, inf*(1+1j)]) 47 with suppress_warnings() as sup: 48 sup.filter(RuntimeWarning, "invalid value encountered in multiply") 49 assert_allclose(spherical_jn(n, x), np.array([0, 0, inf*(1+1j)])) 50 51 def test_spherical_jn_large_arg_1(self): 52 # https://github.com/scipy/scipy/issues/2165 53 # Reference value computed using mpmath, via 54 # besselj(n + mpf(1)/2, z)*sqrt(pi/(2*z)) 55 assert_allclose(spherical_jn(2, 3350.507), -0.00029846226538040747) 56 57 def test_spherical_jn_large_arg_2(self): 58 # https://github.com/scipy/scipy/issues/1641 59 # Reference value computed using mpmath, via 60 # besselj(n + mpf(1)/2, z)*sqrt(pi/(2*z)) 61 assert_allclose(spherical_jn(2, 10000), 3.0590002633029811e-05) 62 63 def test_spherical_jn_at_zero(self): 64 # https://dlmf.nist.gov/10.52.E1 65 # But note that n = 0 is a special case: j0 = sin(x)/x -> 1 66 n = np.array([0, 1, 2, 5, 10, 100]) 67 x = 0 68 assert_allclose(spherical_jn(n, x), np.array([1, 0, 0, 0, 0, 0])) 69 70 71class TestSphericalYn: 72 def test_spherical_yn_exact(self): 73 # https://dlmf.nist.gov/10.49.E5 74 # Note: exact expression is numerically stable only for small 75 # n or z >> n. 76 x = np.array([0.12, 1.23, 12.34, 123.45, 1234.5]) 77 assert_allclose(spherical_yn(2, x), 78 (1/x - 3/x**3)*cos(x) - 3/x**2*sin(x)) 79 80 def test_spherical_yn_recurrence_real(self): 81 # https://dlmf.nist.gov/10.51.E1 82 n = np.array([1, 2, 3, 7, 12]) 83 x = 0.12 84 assert_allclose(spherical_yn(n - 1, x) + spherical_yn(n + 1,x), 85 (2*n + 1)/x*spherical_yn(n, x)) 86 87 def test_spherical_yn_recurrence_complex(self): 88 # https://dlmf.nist.gov/10.51.E1 89 n = np.array([1, 2, 3, 7, 12]) 90 x = 1.1 + 1.5j 91 assert_allclose(spherical_yn(n - 1, x) + spherical_yn(n + 1, x), 92 (2*n + 1)/x*spherical_yn(n, x)) 93 94 def test_spherical_yn_inf_real(self): 95 # https://dlmf.nist.gov/10.52.E3 96 n = 6 97 x = np.array([-inf, inf]) 98 assert_allclose(spherical_yn(n, x), np.array([0, 0])) 99 100 def test_spherical_yn_inf_complex(self): 101 # https://dlmf.nist.gov/10.52.E3 102 n = 7 103 x = np.array([-inf + 0j, inf + 0j, inf*(1+1j)]) 104 with suppress_warnings() as sup: 105 sup.filter(RuntimeWarning, "invalid value encountered in multiply") 106 assert_allclose(spherical_yn(n, x), np.array([0, 0, inf*(1+1j)])) 107 108 def test_spherical_yn_at_zero(self): 109 # https://dlmf.nist.gov/10.52.E2 110 n = np.array([0, 1, 2, 5, 10, 100]) 111 x = 0 112 assert_allclose(spherical_yn(n, x), np.full(n.shape, -inf)) 113 114 def test_spherical_yn_at_zero_complex(self): 115 # Consistently with numpy: 116 # >>> -np.cos(0)/0 117 # -inf 118 # >>> -np.cos(0+0j)/(0+0j) 119 # (-inf + nan*j) 120 n = np.array([0, 1, 2, 5, 10, 100]) 121 x = 0 + 0j 122 assert_allclose(spherical_yn(n, x), np.full(n.shape, nan)) 123 124 125class TestSphericalJnYnCrossProduct: 126 def test_spherical_jn_yn_cross_product_1(self): 127 # https://dlmf.nist.gov/10.50.E3 128 n = np.array([1, 5, 8]) 129 x = np.array([0.1, 1, 10]) 130 left = (spherical_jn(n + 1, x) * spherical_yn(n, x) - 131 spherical_jn(n, x) * spherical_yn(n + 1, x)) 132 right = 1/x**2 133 assert_allclose(left, right) 134 135 def test_spherical_jn_yn_cross_product_2(self): 136 # https://dlmf.nist.gov/10.50.E3 137 n = np.array([1, 5, 8]) 138 x = np.array([0.1, 1, 10]) 139 left = (spherical_jn(n + 2, x) * spherical_yn(n, x) - 140 spherical_jn(n, x) * spherical_yn(n + 2, x)) 141 right = (2*n + 3)/x**3 142 assert_allclose(left, right) 143 144 145class TestSphericalIn: 146 def test_spherical_in_exact(self): 147 # https://dlmf.nist.gov/10.49.E9 148 x = np.array([0.12, 1.23, 12.34, 123.45]) 149 assert_allclose(spherical_in(2, x), 150 (1/x + 3/x**3)*sinh(x) - 3/x**2*cosh(x)) 151 152 def test_spherical_in_recurrence_real(self): 153 # https://dlmf.nist.gov/10.51.E4 154 n = np.array([1, 2, 3, 7, 12]) 155 x = 0.12 156 assert_allclose(spherical_in(n - 1, x) - spherical_in(n + 1,x), 157 (2*n + 1)/x*spherical_in(n, x)) 158 159 def test_spherical_in_recurrence_complex(self): 160 # https://dlmf.nist.gov/10.51.E1 161 n = np.array([1, 2, 3, 7, 12]) 162 x = 1.1 + 1.5j 163 assert_allclose(spherical_in(n - 1, x) - spherical_in(n + 1,x), 164 (2*n + 1)/x*spherical_in(n, x)) 165 166 def test_spherical_in_inf_real(self): 167 # https://dlmf.nist.gov/10.52.E3 168 n = 5 169 x = np.array([-inf, inf]) 170 assert_allclose(spherical_in(n, x), np.array([-inf, inf])) 171 172 def test_spherical_in_inf_complex(self): 173 # https://dlmf.nist.gov/10.52.E5 174 # Ideally, i1n(n, 1j*inf) = 0 and i1n(n, (1+1j)*inf) = (1+1j)*inf, but 175 # this appears impossible to achieve because C99 regards any complex 176 # value with at least one infinite part as a complex infinity, so 177 # 1j*inf cannot be distinguished from (1+1j)*inf. Therefore, nan is 178 # the correct return value. 179 n = 7 180 x = np.array([-inf + 0j, inf + 0j, inf*(1+1j)]) 181 assert_allclose(spherical_in(n, x), np.array([-inf, inf, nan])) 182 183 def test_spherical_in_at_zero(self): 184 # https://dlmf.nist.gov/10.52.E1 185 # But note that n = 0 is a special case: i0 = sinh(x)/x -> 1 186 n = np.array([0, 1, 2, 5, 10, 100]) 187 x = 0 188 assert_allclose(spherical_in(n, x), np.array([1, 0, 0, 0, 0, 0])) 189 190 191class TestSphericalKn: 192 def test_spherical_kn_exact(self): 193 # https://dlmf.nist.gov/10.49.E13 194 x = np.array([0.12, 1.23, 12.34, 123.45]) 195 assert_allclose(spherical_kn(2, x), 196 pi/2*exp(-x)*(1/x + 3/x**2 + 3/x**3)) 197 198 def test_spherical_kn_recurrence_real(self): 199 # https://dlmf.nist.gov/10.51.E4 200 n = np.array([1, 2, 3, 7, 12]) 201 x = 0.12 202 assert_allclose((-1)**(n - 1)*spherical_kn(n - 1, x) - (-1)**(n + 1)*spherical_kn(n + 1,x), 203 (-1)**n*(2*n + 1)/x*spherical_kn(n, x)) 204 205 def test_spherical_kn_recurrence_complex(self): 206 # https://dlmf.nist.gov/10.51.E4 207 n = np.array([1, 2, 3, 7, 12]) 208 x = 1.1 + 1.5j 209 assert_allclose((-1)**(n - 1)*spherical_kn(n - 1, x) - (-1)**(n + 1)*spherical_kn(n + 1,x), 210 (-1)**n*(2*n + 1)/x*spherical_kn(n, x)) 211 212 def test_spherical_kn_inf_real(self): 213 # https://dlmf.nist.gov/10.52.E6 214 n = 5 215 x = np.array([-inf, inf]) 216 assert_allclose(spherical_kn(n, x), np.array([-inf, 0])) 217 218 def test_spherical_kn_inf_complex(self): 219 # https://dlmf.nist.gov/10.52.E6 220 # The behavior at complex infinity depends on the sign of the real 221 # part: if Re(z) >= 0, then the limit is 0; if Re(z) < 0, then it's 222 # z*inf. This distinction cannot be captured, so we return nan. 223 n = 7 224 x = np.array([-inf + 0j, inf + 0j, inf*(1+1j)]) 225 assert_allclose(spherical_kn(n, x), np.array([-inf, 0, nan])) 226 227 def test_spherical_kn_at_zero(self): 228 # https://dlmf.nist.gov/10.52.E2 229 n = np.array([0, 1, 2, 5, 10, 100]) 230 x = 0 231 assert_allclose(spherical_kn(n, x), np.full(n.shape, inf)) 232 233 def test_spherical_kn_at_zero_complex(self): 234 # https://dlmf.nist.gov/10.52.E2 235 n = np.array([0, 1, 2, 5, 10, 100]) 236 x = 0 + 0j 237 assert_allclose(spherical_kn(n, x), np.full(n.shape, nan)) 238 239 240class SphericalDerivativesTestCase: 241 def fundamental_theorem(self, n, a, b): 242 integral, tolerance = quad(lambda z: self.df(n, z), a, b) 243 assert_allclose(integral, 244 self.f(n, b) - self.f(n, a), 245 atol=tolerance) 246 247 @pytest.mark.slow 248 def test_fundamental_theorem_0(self): 249 self.fundamental_theorem(0, 3.0, 15.0) 250 251 @pytest.mark.slow 252 def test_fundamental_theorem_7(self): 253 self.fundamental_theorem(7, 0.5, 1.2) 254 255 256class TestSphericalJnDerivatives(SphericalDerivativesTestCase): 257 def f(self, n, z): 258 return spherical_jn(n, z) 259 260 def df(self, n, z): 261 return spherical_jn(n, z, derivative=True) 262 263 def test_spherical_jn_d_zero(self): 264 n = np.array([0, 1, 2, 3, 7, 15]) 265 assert_allclose(spherical_jn(n, 0, derivative=True), 266 np.array([0, 1/3, 0, 0, 0, 0])) 267 268 269class TestSphericalYnDerivatives(SphericalDerivativesTestCase): 270 def f(self, n, z): 271 return spherical_yn(n, z) 272 273 def df(self, n, z): 274 return spherical_yn(n, z, derivative=True) 275 276 277class TestSphericalInDerivatives(SphericalDerivativesTestCase): 278 def f(self, n, z): 279 return spherical_in(n, z) 280 281 def df(self, n, z): 282 return spherical_in(n, z, derivative=True) 283 284 def test_spherical_in_d_zero(self): 285 n = np.array([1, 2, 3, 7, 15]) 286 assert_allclose(spherical_in(n, 0, derivative=True), 287 np.zeros(5)) 288 289 290class TestSphericalKnDerivatives(SphericalDerivativesTestCase): 291 def f(self, n, z): 292 return spherical_kn(n, z) 293 294 def df(self, n, z): 295 return spherical_kn(n, z, derivative=True) 296 297 298class TestSphericalOld: 299 # These are tests from the TestSpherical class of test_basic.py, 300 # rewritten to use spherical_* instead of sph_* but otherwise unchanged. 301 302 def test_sph_in(self): 303 # This test reproduces test_basic.TestSpherical.test_sph_in. 304 i1n = np.empty((2,2)) 305 x = 0.2 306 307 i1n[0][0] = spherical_in(0, x) 308 i1n[0][1] = spherical_in(1, x) 309 i1n[1][0] = spherical_in(0, x, derivative=True) 310 i1n[1][1] = spherical_in(1, x, derivative=True) 311 312 inp0 = (i1n[0][1]) 313 inp1 = (i1n[0][0] - 2.0/0.2 * i1n[0][1]) 314 assert_array_almost_equal(i1n[0],np.array([1.0066800127054699381, 315 0.066933714568029540839]),12) 316 assert_array_almost_equal(i1n[1],[inp0,inp1],12) 317 318 def test_sph_in_kn_order0(self): 319 x = 1. 320 sph_i0 = np.empty((2,)) 321 sph_i0[0] = spherical_in(0, x) 322 sph_i0[1] = spherical_in(0, x, derivative=True) 323 sph_i0_expected = np.array([np.sinh(x)/x, 324 np.cosh(x)/x-np.sinh(x)/x**2]) 325 assert_array_almost_equal(r_[sph_i0], sph_i0_expected) 326 327 sph_k0 = np.empty((2,)) 328 sph_k0[0] = spherical_kn(0, x) 329 sph_k0[1] = spherical_kn(0, x, derivative=True) 330 sph_k0_expected = np.array([0.5*pi*exp(-x)/x, 331 -0.5*pi*exp(-x)*(1/x+1/x**2)]) 332 assert_array_almost_equal(r_[sph_k0], sph_k0_expected) 333 334 def test_sph_jn(self): 335 s1 = np.empty((2,3)) 336 x = 0.2 337 338 s1[0][0] = spherical_jn(0, x) 339 s1[0][1] = spherical_jn(1, x) 340 s1[0][2] = spherical_jn(2, x) 341 s1[1][0] = spherical_jn(0, x, derivative=True) 342 s1[1][1] = spherical_jn(1, x, derivative=True) 343 s1[1][2] = spherical_jn(2, x, derivative=True) 344 345 s10 = -s1[0][1] 346 s11 = s1[0][0]-2.0/0.2*s1[0][1] 347 s12 = s1[0][1]-3.0/0.2*s1[0][2] 348 assert_array_almost_equal(s1[0],[0.99334665397530607731, 349 0.066400380670322230863, 350 0.0026590560795273856680],12) 351 assert_array_almost_equal(s1[1],[s10,s11,s12],12) 352 353 def test_sph_kn(self): 354 kn = np.empty((2,3)) 355 x = 0.2 356 357 kn[0][0] = spherical_kn(0, x) 358 kn[0][1] = spherical_kn(1, x) 359 kn[0][2] = spherical_kn(2, x) 360 kn[1][0] = spherical_kn(0, x, derivative=True) 361 kn[1][1] = spherical_kn(1, x, derivative=True) 362 kn[1][2] = spherical_kn(2, x, derivative=True) 363 364 kn0 = -kn[0][1] 365 kn1 = -kn[0][0]-2.0/0.2*kn[0][1] 366 kn2 = -kn[0][1]-3.0/0.2*kn[0][2] 367 assert_array_almost_equal(kn[0],[6.4302962978445670140, 368 38.581777787067402086, 369 585.15696310385559829],12) 370 assert_array_almost_equal(kn[1],[kn0,kn1,kn2],9) 371 372 def test_sph_yn(self): 373 sy1 = spherical_yn(2, 0.2) 374 sy2 = spherical_yn(0, 0.2) 375 assert_almost_equal(sy1,-377.52483,5) # previous values in the system 376 assert_almost_equal(sy2,-4.9003329,5) 377 sphpy = (spherical_yn(0, 0.2) - 2*spherical_yn(2, 0.2))/3 378 sy3 = spherical_yn(1, 0.2, derivative=True) 379 assert_almost_equal(sy3,sphpy,4) # compare correct derivative val. (correct =-system val). 380