1import numpy as np 2 3from ._ufuncs import _ellip_harm 4from ._ellip_harm_2 import _ellipsoid, _ellipsoid_norm 5 6 7def ellip_harm(h2, k2, n, p, s, signm=1, signn=1): 8 r""" 9 Ellipsoidal harmonic functions E^p_n(l) 10 11 These are also known as Lame functions of the first kind, and are 12 solutions to the Lame equation: 13 14 .. math:: (s^2 - h^2)(s^2 - k^2)E''(s) + s(2s^2 - h^2 - k^2)E'(s) + (a - q s^2)E(s) = 0 15 16 where :math:`q = (n+1)n` and :math:`a` is the eigenvalue (not 17 returned) corresponding to the solutions. 18 19 Parameters 20 ---------- 21 h2 : float 22 ``h**2`` 23 k2 : float 24 ``k**2``; should be larger than ``h**2`` 25 n : int 26 Degree 27 s : float 28 Coordinate 29 p : int 30 Order, can range between [1,2n+1] 31 signm : {1, -1}, optional 32 Sign of prefactor of functions. Can be +/-1. See Notes. 33 signn : {1, -1}, optional 34 Sign of prefactor of functions. Can be +/-1. See Notes. 35 36 Returns 37 ------- 38 E : float 39 the harmonic :math:`E^p_n(s)` 40 41 See Also 42 -------- 43 ellip_harm_2, ellip_normal 44 45 Notes 46 ----- 47 The geometric interpretation of the ellipsoidal functions is 48 explained in [2]_, [3]_, [4]_. The `signm` and `signn` arguments control the 49 sign of prefactors for functions according to their type:: 50 51 K : +1 52 L : signm 53 M : signn 54 N : signm*signn 55 56 .. versionadded:: 0.15.0 57 58 References 59 ---------- 60 .. [1] Digital Library of Mathematical Functions 29.12 61 https://dlmf.nist.gov/29.12 62 .. [2] Bardhan and Knepley, "Computational science and 63 re-discovery: open-source implementations of 64 ellipsoidal harmonics for problems in potential theory", 65 Comput. Sci. Disc. 5, 014006 (2012) 66 :doi:`10.1088/1749-4699/5/1/014006`. 67 .. [3] David J.and Dechambre P, "Computation of Ellipsoidal 68 Gravity Field Harmonics for small solar system bodies" 69 pp. 30-36, 2000 70 .. [4] George Dassios, "Ellipsoidal Harmonics: Theory and Applications" 71 pp. 418, 2012 72 73 Examples 74 -------- 75 >>> from scipy.special import ellip_harm 76 >>> w = ellip_harm(5,8,1,1,2.5) 77 >>> w 78 2.5 79 80 Check that the functions indeed are solutions to the Lame equation: 81 82 >>> from scipy.interpolate import UnivariateSpline 83 >>> def eigenvalue(f, df, ddf): 84 ... r = ((s**2 - h**2)*(s**2 - k**2)*ddf + s*(2*s**2 - h**2 - k**2)*df - n*(n+1)*s**2*f)/f 85 ... return -r.mean(), r.std() 86 >>> s = np.linspace(0.1, 10, 200) 87 >>> k, h, n, p = 8.0, 2.2, 3, 2 88 >>> E = ellip_harm(h**2, k**2, n, p, s) 89 >>> E_spl = UnivariateSpline(s, E) 90 >>> a, a_err = eigenvalue(E_spl(s), E_spl(s,1), E_spl(s,2)) 91 >>> a, a_err 92 (583.44366156701483, 6.4580890640310646e-11) 93 94 """ 95 return _ellip_harm(h2, k2, n, p, s, signm, signn) 96 97 98_ellip_harm_2_vec = np.vectorize(_ellipsoid, otypes='d') 99 100 101def ellip_harm_2(h2, k2, n, p, s): 102 r""" 103 Ellipsoidal harmonic functions F^p_n(l) 104 105 These are also known as Lame functions of the second kind, and are 106 solutions to the Lame equation: 107 108 .. math:: (s^2 - h^2)(s^2 - k^2)F''(s) + s(2s^2 - h^2 - k^2)F'(s) + (a - q s^2)F(s) = 0 109 110 where :math:`q = (n+1)n` and :math:`a` is the eigenvalue (not 111 returned) corresponding to the solutions. 112 113 Parameters 114 ---------- 115 h2 : float 116 ``h**2`` 117 k2 : float 118 ``k**2``; should be larger than ``h**2`` 119 n : int 120 Degree. 121 p : int 122 Order, can range between [1,2n+1]. 123 s : float 124 Coordinate 125 126 Returns 127 ------- 128 F : float 129 The harmonic :math:`F^p_n(s)` 130 131 Notes 132 ----- 133 Lame functions of the second kind are related to the functions of the first kind: 134 135 .. math:: 136 137 F^p_n(s)=(2n + 1)E^p_n(s)\int_{0}^{1/s}\frac{du}{(E^p_n(1/u))^2\sqrt{(1-u^2k^2)(1-u^2h^2)}} 138 139 .. versionadded:: 0.15.0 140 141 See Also 142 -------- 143 ellip_harm, ellip_normal 144 145 Examples 146 -------- 147 >>> from scipy.special import ellip_harm_2 148 >>> w = ellip_harm_2(5,8,2,1,10) 149 >>> w 150 0.00108056853382 151 152 """ 153 with np.errstate(all='ignore'): 154 return _ellip_harm_2_vec(h2, k2, n, p, s) 155 156 157def _ellip_normal_vec(h2, k2, n, p): 158 return _ellipsoid_norm(h2, k2, n, p) 159 160 161_ellip_normal_vec = np.vectorize(_ellip_normal_vec, otypes='d') 162 163 164def ellip_normal(h2, k2, n, p): 165 r""" 166 Ellipsoidal harmonic normalization constants gamma^p_n 167 168 The normalization constant is defined as 169 170 .. math:: 171 172 \gamma^p_n=8\int_{0}^{h}dx\int_{h}^{k}dy\frac{(y^2-x^2)(E^p_n(y)E^p_n(x))^2}{\sqrt((k^2-y^2)(y^2-h^2)(h^2-x^2)(k^2-x^2)} 173 174 Parameters 175 ---------- 176 h2 : float 177 ``h**2`` 178 k2 : float 179 ``k**2``; should be larger than ``h**2`` 180 n : int 181 Degree. 182 p : int 183 Order, can range between [1,2n+1]. 184 185 Returns 186 ------- 187 gamma : float 188 The normalization constant :math:`\gamma^p_n` 189 190 See Also 191 -------- 192 ellip_harm, ellip_harm_2 193 194 Notes 195 ----- 196 .. versionadded:: 0.15.0 197 198 Examples 199 -------- 200 >>> from scipy.special import ellip_normal 201 >>> w = ellip_normal(5,8,3,7) 202 >>> w 203 1723.38796997 204 205 """ 206 with np.errstate(all='ignore'): 207 return _ellip_normal_vec(h2, k2, n, p) 208