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