1#
2# Tests for the Ellipsoidal Harmonic Function,
3# Distributed under the same license as SciPy itself.
4#
5
6import numpy as np
7from numpy.testing import (assert_equal, assert_almost_equal, assert_allclose,
8                           assert_, suppress_warnings)
9from scipy.special._testutils import assert_func_equal
10from scipy.special import ellip_harm, ellip_harm_2, ellip_normal
11from scipy.integrate import IntegrationWarning
12from numpy import sqrt, pi
13
14
15def test_ellip_potential():
16    def change_coefficient(lambda1, mu, nu, h2, k2):
17        x = sqrt(lambda1**2*mu**2*nu**2/(h2*k2))
18        y = sqrt((lambda1**2 - h2)*(mu**2 - h2)*(h2 - nu**2)/(h2*(k2 - h2)))
19        z = sqrt((lambda1**2 - k2)*(k2 - mu**2)*(k2 - nu**2)/(k2*(k2 - h2)))
20        return x, y, z
21
22    def solid_int_ellip(lambda1, mu, nu, n, p, h2, k2):
23        return (ellip_harm(h2, k2, n, p, lambda1)*ellip_harm(h2, k2, n, p, mu)
24               * ellip_harm(h2, k2, n, p, nu))
25
26    def solid_int_ellip2(lambda1, mu, nu, n, p, h2, k2):
27        return (ellip_harm_2(h2, k2, n, p, lambda1)
28                * ellip_harm(h2, k2, n, p, mu)*ellip_harm(h2, k2, n, p, nu))
29
30    def summation(lambda1, mu1, nu1, lambda2, mu2, nu2, h2, k2):
31        tol = 1e-8
32        sum1 = 0
33        for n in range(20):
34            xsum = 0
35            for p in range(1, 2*n+2):
36                xsum += (4*pi*(solid_int_ellip(lambda2, mu2, nu2, n, p, h2, k2)
37                    * solid_int_ellip2(lambda1, mu1, nu1, n, p, h2, k2)) /
38                    (ellip_normal(h2, k2, n, p)*(2*n + 1)))
39            if abs(xsum) < 0.1*tol*abs(sum1):
40                break
41            sum1 += xsum
42        return sum1, xsum
43
44    def potential(lambda1, mu1, nu1, lambda2, mu2, nu2, h2, k2):
45        x1, y1, z1 = change_coefficient(lambda1, mu1, nu1, h2, k2)
46        x2, y2, z2 = change_coefficient(lambda2, mu2, nu2, h2, k2)
47        res = sqrt((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)
48        return 1/res
49
50    pts = [
51        (120, sqrt(19), 2, 41, sqrt(17), 2, 15, 25),
52        (120, sqrt(16), 3.2, 21, sqrt(11), 2.9, 11, 20),
53       ]
54
55    with suppress_warnings() as sup:
56        sup.filter(IntegrationWarning, "The occurrence of roundoff error")
57        sup.filter(IntegrationWarning, "The maximum number of subdivisions")
58
59        for p in pts:
60            err_msg = repr(p)
61            exact = potential(*p)
62            result, last_term = summation(*p)
63            assert_allclose(exact, result, atol=0, rtol=1e-8, err_msg=err_msg)
64            assert_(abs(result - exact) < 10*abs(last_term), err_msg)
65
66
67def test_ellip_norm():
68
69    def G01(h2, k2):
70        return 4*pi
71
72    def G11(h2, k2):
73        return 4*pi*h2*k2/3
74
75    def G12(h2, k2):
76        return 4*pi*h2*(k2 - h2)/3
77
78    def G13(h2, k2):
79        return 4*pi*k2*(k2 - h2)/3
80
81    def G22(h2, k2):
82        res = (2*(h2**4 + k2**4) - 4*h2*k2*(h2**2 + k2**2) + 6*h2**2*k2**2 +
83        sqrt(h2**2 + k2**2 - h2*k2)*(-2*(h2**3 + k2**3) + 3*h2*k2*(h2 + k2)))
84        return 16*pi/405*res
85
86    def G21(h2, k2):
87        res = (2*(h2**4 + k2**4) - 4*h2*k2*(h2**2 + k2**2) + 6*h2**2*k2**2
88        + sqrt(h2**2 + k2**2 - h2*k2)*(2*(h2**3 + k2**3) - 3*h2*k2*(h2 + k2)))
89        return 16*pi/405*res
90
91    def G23(h2, k2):
92        return 4*pi*h2**2*k2*(k2 - h2)/15
93
94    def G24(h2, k2):
95        return 4*pi*h2*k2**2*(k2 - h2)/15
96
97    def G25(h2, k2):
98        return 4*pi*h2*k2*(k2 - h2)**2/15
99
100    def G32(h2, k2):
101        res = (16*(h2**4 + k2**4) - 36*h2*k2*(h2**2 + k2**2) + 46*h2**2*k2**2
102        + sqrt(4*(h2**2 + k2**2) - 7*h2*k2)*(-8*(h2**3 + k2**3) +
103        11*h2*k2*(h2 + k2)))
104        return 16*pi/13125*k2*h2*res
105
106    def G31(h2, k2):
107        res = (16*(h2**4 + k2**4) - 36*h2*k2*(h2**2 + k2**2) + 46*h2**2*k2**2
108        + sqrt(4*(h2**2 + k2**2) - 7*h2*k2)*(8*(h2**3 + k2**3) -
109        11*h2*k2*(h2 + k2)))
110        return 16*pi/13125*h2*k2*res
111
112    def G34(h2, k2):
113        res = (6*h2**4 + 16*k2**4 - 12*h2**3*k2 - 28*h2*k2**3 + 34*h2**2*k2**2
114        + sqrt(h2**2 + 4*k2**2 - h2*k2)*(-6*h2**3 - 8*k2**3 + 9*h2**2*k2 +
115                                            13*h2*k2**2))
116        return 16*pi/13125*h2*(k2 - h2)*res
117
118    def G33(h2, k2):
119        res = (6*h2**4 + 16*k2**4 - 12*h2**3*k2 - 28*h2*k2**3 + 34*h2**2*k2**2
120        + sqrt(h2**2 + 4*k2**2 - h2*k2)*(6*h2**3 + 8*k2**3 - 9*h2**2*k2 -
121        13*h2*k2**2))
122        return 16*pi/13125*h2*(k2 - h2)*res
123
124    def G36(h2, k2):
125        res = (16*h2**4 + 6*k2**4 - 28*h2**3*k2 - 12*h2*k2**3 + 34*h2**2*k2**2
126        + sqrt(4*h2**2 + k2**2 - h2*k2)*(-8*h2**3 - 6*k2**3 + 13*h2**2*k2 +
127        9*h2*k2**2))
128        return 16*pi/13125*k2*(k2 - h2)*res
129
130    def G35(h2, k2):
131        res = (16*h2**4 + 6*k2**4 - 28*h2**3*k2 - 12*h2*k2**3 + 34*h2**2*k2**2
132        + sqrt(4*h2**2 + k2**2 - h2*k2)*(8*h2**3 + 6*k2**3 - 13*h2**2*k2 -
133        9*h2*k2**2))
134        return 16*pi/13125*k2*(k2 - h2)*res
135
136    def G37(h2, k2):
137        return 4*pi*h2**2*k2**2*(k2 - h2)**2/105
138
139    known_funcs = {(0, 1): G01, (1, 1): G11, (1, 2): G12, (1, 3): G13,
140                   (2, 1): G21, (2, 2): G22, (2, 3): G23, (2, 4): G24,
141                   (2, 5): G25, (3, 1): G31, (3, 2): G32, (3, 3): G33,
142                   (3, 4): G34, (3, 5): G35, (3, 6): G36, (3, 7): G37}
143
144    def _ellip_norm(n, p, h2, k2):
145        func = known_funcs[n, p]
146        return func(h2, k2)
147    _ellip_norm = np.vectorize(_ellip_norm)
148
149    def ellip_normal_known(h2, k2, n, p):
150        return _ellip_norm(n, p, h2, k2)
151
152    # generate both large and small h2 < k2 pairs
153    np.random.seed(1234)
154    h2 = np.random.pareto(0.5, size=1)
155    k2 = h2 * (1 + np.random.pareto(0.5, size=h2.size))
156
157    points = []
158    for n in range(4):
159        for p in range(1, 2*n+2):
160            points.append((h2, k2, np.full(h2.size, n), np.full(h2.size, p)))
161    points = np.array(points)
162    with suppress_warnings() as sup:
163        sup.filter(IntegrationWarning, "The occurrence of roundoff error")
164        assert_func_equal(ellip_normal, ellip_normal_known, points, rtol=1e-12)
165
166
167def test_ellip_harm_2():
168
169    def I1(h2, k2, s):
170        res = (ellip_harm_2(h2, k2, 1, 1, s)/(3 * ellip_harm(h2, k2, 1, 1, s))
171        + ellip_harm_2(h2, k2, 1, 2, s)/(3 * ellip_harm(h2, k2, 1, 2, s)) +
172        ellip_harm_2(h2, k2, 1, 3, s)/(3 * ellip_harm(h2, k2, 1, 3, s)))
173        return res
174
175    with suppress_warnings() as sup:
176        sup.filter(IntegrationWarning, "The occurrence of roundoff error")
177        assert_almost_equal(I1(5, 8, 10), 1/(10*sqrt((100-5)*(100-8))))
178
179        # Values produced by code from arXiv:1204.0267
180        assert_almost_equal(ellip_harm_2(5, 8, 2, 1, 10), 0.00108056853382)
181        assert_almost_equal(ellip_harm_2(5, 8, 2, 2, 10), 0.00105820513809)
182        assert_almost_equal(ellip_harm_2(5, 8, 2, 3, 10), 0.00106058384743)
183        assert_almost_equal(ellip_harm_2(5, 8, 2, 4, 10), 0.00106774492306)
184        assert_almost_equal(ellip_harm_2(5, 8, 2, 5, 10), 0.00107976356454)
185
186
187def test_ellip_harm():
188
189    def E01(h2, k2, s):
190        return 1
191
192    def E11(h2, k2, s):
193        return s
194
195    def E12(h2, k2, s):
196        return sqrt(abs(s*s - h2))
197
198    def E13(h2, k2, s):
199        return sqrt(abs(s*s - k2))
200
201    def E21(h2, k2, s):
202        return s*s - 1/3*((h2 + k2) + sqrt(abs((h2 + k2)*(h2 + k2)-3*h2*k2)))
203
204    def E22(h2, k2, s):
205        return s*s - 1/3*((h2 + k2) - sqrt(abs((h2 + k2)*(h2 + k2)-3*h2*k2)))
206
207    def E23(h2, k2, s):
208        return s * sqrt(abs(s*s - h2))
209
210    def E24(h2, k2, s):
211        return s * sqrt(abs(s*s - k2))
212
213    def E25(h2, k2, s):
214        return sqrt(abs((s*s - h2)*(s*s - k2)))
215
216    def E31(h2, k2, s):
217        return s*s*s - (s/5)*(2*(h2 + k2) + sqrt(4*(h2 + k2)*(h2 + k2) -
218        15*h2*k2))
219
220    def E32(h2, k2, s):
221        return s*s*s - (s/5)*(2*(h2 + k2) - sqrt(4*(h2 + k2)*(h2 + k2) -
222        15*h2*k2))
223
224    def E33(h2, k2, s):
225        return sqrt(abs(s*s - h2))*(s*s - 1/5*((h2 + 2*k2) + sqrt(abs((h2 +
226        2*k2)*(h2 + 2*k2) - 5*h2*k2))))
227
228    def E34(h2, k2, s):
229        return sqrt(abs(s*s - h2))*(s*s - 1/5*((h2 + 2*k2) - sqrt(abs((h2 +
230        2*k2)*(h2 + 2*k2) - 5*h2*k2))))
231
232    def E35(h2, k2, s):
233        return sqrt(abs(s*s - k2))*(s*s - 1/5*((2*h2 + k2) + sqrt(abs((2*h2
234        + k2)*(2*h2 + k2) - 5*h2*k2))))
235
236    def E36(h2, k2, s):
237        return sqrt(abs(s*s - k2))*(s*s - 1/5*((2*h2 + k2) - sqrt(abs((2*h2
238        + k2)*(2*h2 + k2) - 5*h2*k2))))
239
240    def E37(h2, k2, s):
241        return s * sqrt(abs((s*s - h2)*(s*s - k2)))
242
243    assert_equal(ellip_harm(5, 8, 1, 2, 2.5, 1, 1),
244    ellip_harm(5, 8, 1, 2, 2.5))
245
246    known_funcs = {(0, 1): E01, (1, 1): E11, (1, 2): E12, (1, 3): E13,
247                   (2, 1): E21, (2, 2): E22, (2, 3): E23, (2, 4): E24,
248                   (2, 5): E25, (3, 1): E31, (3, 2): E32, (3, 3): E33,
249                   (3, 4): E34, (3, 5): E35, (3, 6): E36, (3, 7): E37}
250
251    point_ref = []
252
253    def ellip_harm_known(h2, k2, n, p, s):
254        for i in range(h2.size):
255            func = known_funcs[(int(n[i]), int(p[i]))]
256            point_ref.append(func(h2[i], k2[i], s[i]))
257        return point_ref
258
259    np.random.seed(1234)
260    h2 = np.random.pareto(0.5, size=30)
261    k2 = h2*(1 + np.random.pareto(0.5, size=h2.size))
262    s = np.random.pareto(0.5, size=h2.size)
263    points = []
264    for i in range(h2.size):
265        for n in range(4):
266            for p in range(1, 2*n+2):
267                points.append((h2[i], k2[i], n, p, s[i]))
268    points = np.array(points)
269    assert_func_equal(ellip_harm, ellip_harm_known, points, rtol=1e-12)
270
271
272def test_ellip_harm_invalid_p():
273    # Regression test. This should return nan.
274    n = 4
275    # Make p > 2*n + 1.
276    p = 2*n + 2
277    result = ellip_harm(0.5, 2.0, n, p, 0.2)
278    assert np.isnan(result)
279