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