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