Coverage for pygeodesy/elliptic.py : 95%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
C++ class U{EllipticFunction <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1EllipticFunction.html>} into pure Python class L{Elliptic}.
Python method names follow the C++ member functions, except:
- member functions I{without arguments} are mapped to Python properties prefixed with C{"c"}, for example C{E()} is property C{cE},
- member functions with 1 or 3 arguments are renamed to Python methods starting with an C{"f"}, example C{E(psi)} to C{fE(psi)} and C{E(sn, cn, dn)} to C{fE(sn, cn, dn)},
- other Python method names conventionally start with a lower-case letter or an underscore if private.
Following is a copy of Karney's U{EllipticFunction.hpp <https://GeographicLib.SourceForge.io/html/EllipticFunction_8hpp_source.html>} file C{Header}.
Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2017) and licensed under the MIT/X11 License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io/>} documentation.
B{Elliptic integrals and functions.}
This provides the elliptic functions and integrals needed for C{Ellipsoid}, C{GeodesicExact}, and C{TransverseMercatorExact}. Two categories of function are provided:
- functions to compute U{symmetric elliptic integrals <https://DLMF.NIST.gov/19.16.i>}
- methods to compute U{Legrendre's elliptic integrals <https://DLMF.NIST.gov/19.2.ii>} and U{Jacobi elliptic functions<https://DLMF.NIST.gov/22.2>}.
In the latter case, an object is constructed giving the modulus C{k} (and optionally the parameter C{alpha}). The modulus (and parameter) are always passed as squares which allows C{k} to be pure imaginary. (Confusingly, Abramowitz and Stegun call C{m = k**2} the "parameter" and C{n = alpha**2} the "characteristic".)
In geodesic applications, it is convenient to separate the incomplete integrals into secular and periodic components, e.g.
I{C{E(phi, k) = (2 E(k) / pi) [ phi + delta E(phi, k) ]}}
where I{C{delta E(phi, k)}} is an odd periodic function with period I{C{pi}}.
The computation of the elliptic integrals uses the algorithms given in U{B. C. Carlson, Computation of real or complex elliptic integrals <https://DOI.org/10.1007/BF02198293>}, Numerical Algorithms 10, 13--26 (1995) with the additional optimizations given U{here <https://DLMF.NIST.gov/19.36.i>}.
The computation of the Jacobi elliptic functions uses the algorithm given in U{R. Bulirsch, Numerical Calculation of Elliptic Integrals and Elliptic Functions <https://DOI.org/10.1007/BF01397975>}, Numerische Mathematik 7, 78--90 (1965).
The notation follows U{NIST Digital Library of Mathematical Functions <https://DLMF.NIST.gov>} chapters U{19<https://DLMF.NIST.gov/19>} and U{22<https://DLMF.NIST.gov/22>}. ''' # make sure int/int division yields float quotient raise ImportError('%s 1/2 == %d' % ('division', division))
INF, map2 sincos2, sincos2d
floor, sin, sqrt, tanh
'''Elliptic integral, function, convergence or other L{Elliptic} issue. '''
'''Elliptic integrals and functions.
@see: Karney's U{Detailed Description<https://GeographicLib.SourceForge.io/ html/classGeographicLib_1_1EllipticFunction.html#details>}. '''
'''Constructor, specifying the C{modulus} and C{parameter}.
@keyword k2: Modulus squared (C{scalar}, 0 <= k^2 <= 1). @keyword alpha2: Parameter squared (C{scalar}, 0 <= α^2 <= 1). @keyword kp2: Complementary modulus squared (C{float}, k'^2 >= 0). @keyword alphap2: Complementary parameter squared (C{float}, α'^2 >= 0).
@see: Method L{reset} for further details.
@note: If only elliptic integrals of the first and second kinds are needed, then set B{C{alpha2}} = 0 (the default value). In that case, we have Π(φ, 0, k) = F(φ, k), G(φ, 0, k) = E(φ, k), and H(φ, 0, k) = F(φ, k) - D(φ, k). '''
def alpha2(self): '''Get α^2, the square of the parameter (C{float}). '''
def alphap2(self): '''Get α'^2, the square of the complementary parameter (C{float}). '''
def cD(self): '''Get Jahnke's complete integral C{D(k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E6>}. '''
def cE(self): '''Get the complete integral of the second kind C{E(k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}. '''
def cG(self): '''Get Legendre's complete geodesic longitude integral C{G(alpha2, k)} (C{float}). '''
def cH(self): '''Get Cayley's complete geodesic longitude difference integral C{H(alpha2, k)} (C{float}). '''
def cK(self): '''Get the complete integral of the first kind C{K(k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}. '''
def cKE(self): '''Get the difference between the complete integrals of the first and second kinds, C{K(k) − E(k)} (C{float}). '''
def cPi(self): '''Get the complete integral of the third kind C{Pi(alpha2, k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}. '''
'''The periodic Jahnke's incomplete elliptic integral.
@param sn: sin(φ). @param cn: cos(φ). @param dn: sqrt(1 − k2 sin(2φ)).
@return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}).
@raise EllipticError: Invalid invokation or no convergence. '''
'''The periodic incomplete integral of the second kind.
@param sn: sin(φ). @param cn: cos(φ). @param dn: sqrt(1 − k2 sin(2φ)).
@return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}).
@raise EllipticError: Invalid invokation or no convergence. '''
'''The periodic inverse of the incomplete integral of the second kind.
@param stau: sin(τ) @param ctau: cos(τ)
@return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}).
@raise EllipticError: No convergence. ''' # Function is periodic with period pi
'''The periodic incomplete integral of the first kind.
@param sn: sin(φ). @param cn: cos(φ). @param dn: sqrt(1 − k2 sin(2φ)).
@return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}).
@raise EllipticError: Invalid invokation or no convergence. '''
'''Legendre's periodic geodesic longitude integral.
@param sn: sin(φ). @param cn: cos(φ). @param dn: sqrt(1 − k2 sin(2φ)).
@return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}).
@raise EllipticError: Invalid invokation or no convergence. '''
'''Cayley's periodic geodesic longitude difference integral.
@param sn: sin(φ). @param cn: cos(φ). @param dn: sqrt(1 − k2 sin(2φ)).
@return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}).
@raise EllipticError: Invalid invokation or no convergence. '''
'''The periodic incomplete integral of the third kind.
@param sn: sin(φ). @param cn: cos(φ). @param dn: sqrt(1 − k2 sin(2φ)).
@return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ (C{float}).
@raise EllipticError: Invalid invokation or no convergence. '''
'''(INTERNAL) Helper for C{.deltaD} thru C{.deltaPi}. ''' raise EllipticError('%s invalid: %s%r' % ('invokation', t, (sn, cn, dn)))
cn, sn = -cn, -sn
def eps(self): '''Get epsilon (C{float}). '''
'''Jahnke's incomplete elliptic integral in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sin(φ). @keyword cn: C{None} or cos(φ). @keyword dn: C{None} or sqrt(1 − k2 sin(2φ)).
@return: D(φ, k) as though φ ∈ (−π, π] (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}.
@raise EllipticError: Invalid invokation or no convergence. '''
self.deltaD, _fD)
'''The C{Delta} amplitude function.
@param sn: sin(φ). @param cn: cos(φ).
@return: C{sqrt(1 − k2 C{sin}(2φ))} (C{float}). ''' (k2 * cn**2 + self.kp2))
'''The incomplete integral of the second kind in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sin(φ). @keyword cn: C{None} or cos(φ). @keyword dn: C{None} or sqrt(1 − k2 sin(2φ)).
@return: E(φ, k) as though φ ∈ (−π, π] (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}.
@raise EllipticError: Invalid invokation or no convergence. ''' # Carlson, eq. 4.6 and <https://DLMF.NIST.gov/19.25.E9> ei = _RF(cn2, dn2, 1) - k2 * sn2 * _RD_3(cn2, dn2, 1) # <https://DLMF.NIST.gov/19.25.E10> kp2 * k2 * sn2 * _RD_3(cn2, 1, dn2), k2 * abs(cn) / dn) else: # <https://DLMF.NIST.gov/19.25.E11> ei = dn / abs(cn) - kp2 * sn2 * _RD_3(dn2, 1, cn2)
self.deltaE, _fE)
'''The incomplete integral of the second kind with the argument given in degrees.
@param deg: Angle (C{degrees}).
@return: E(π B{C{deg}}/180, k) (C{float}).
@raise EllipticError: No convergence. '''
'''The inverse of the incomplete integral of the second kind.
@param x: Argument (C{float}).
@return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}} (C{float}).
@raise EllipticError: No convergence. ''' # linear approximation # first order correction # For kp2 close to zero use asin(x/.cE) or J. P. Boyd, # Applied Math. and Computation 218, 7005-7013 (2012) # <https://DOI.org/10.1016/j.amc.2011.12.021> raise EllipticError('no %s%r convergence' % ('fEinv', (x,)))
'''The incomplete integral of the first kind in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sin(φ). @keyword cn: C{None} or cos(φ). @keyword dn: C{None} or sqrt(1 − k2 sin(2φ)).
@return: F(φ, k) as though φ ∈ (−π, π] (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}.
@raise EllipticError: Invalid invokation or no convergence. '''
self.deltaF, _fF)
'''Legendre's geodesic longitude integral in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sin(φ). @keyword cn: C{None} or cos(φ). @keyword dn: C{None} or sqrt(1 − k2 sin(2φ)).
@return: G(φ, k) as though φ ∈ (−π, π] (C{float}).
@raise EllipticError: Invalid invokation or no convergence.
@note: Legendre expresses the longitude of a point on the geodesic in terms of this combination of elliptic integrals in U{Exercices de Calcul Intégral, Vol 1 (1811), p 181<https://Books. Google.com/books?id=riIOAAAAQAAJ&pg=PA181>}.
@see: U{Geodesics in terms of elliptic integrals<https:// GeographicLib.SourceForge.io/html/geodesic.html#geodellip>} for the expression for the longitude in terms of this function. ''' self.cG, self.deltaG)
'''Cayley's geodesic longitude difference integral in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sin(φ). @keyword cn: C{None} or cos(φ). @keyword dn: C{None} or sqrt(1 − k2 sin(2φ)).
@return: H(φ, k) as though φ ∈ (−π, π] (C{float}).
@raise EllipticError: Invalid invokation or no convergence.
@note: Cayley expresses the longitude difference of a point on the geodesic in terms of this combination of elliptic integrals in U{Phil. Mag. B{40} (1870), p 333 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}.
@see: U{Geodesics in terms of elliptic integrals<https:// GeographicLib.SourceForge.io/html/geodesic.html#geodellip>} for the expression for the longitude in terms of this function. ''' self.cH, self.deltaH)
'''The incomplete integral of the third kind in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sin(φ). @keyword cn: C{None} or cos(φ). @keyword dn: C{None} or sqrt(1 − k2 sin(2φ)).
@return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}).
@raise EllipticError: Invalid invokation or no convergence. ''' self.cPi, self.deltaPi)
'''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}. ''' _RJ_3(cn2, dn2, 1, cn2 + self.alphap2 * sn2))
'''(INTERNAL) Helper for C{f.D}, C{.fE}, C{.fF} and C{._fXa}. ''' return (deltaX(sn, cn, dn) + phi) * cX / PI_2 # fall through raise EllipticError('%s invalid: %s%r' % ('invokation', t, (sn, cn, dn)))
elif cn < 0: xi = 2 * cX - fX(sn, cn, dn) else: xi = cX
def k2(self): '''Get k^2, the square of the modulus (C{float}). '''
def kp2(self): '''Get k'^2, the square of the complementary modulus (C{float}). '''
'''Reset the modulus and parameter.
@keyword k2: modulus squared (C{float}, 0 <= k^2<= 1). @keyword alpha2: parameter (C{float}, 0 <= α^2 <= 1). @keyword kp2: complementary modulus squared (C{float}, k'^2 >= 0). @keyword alphap2: complementary parameter squared (C{float}, α'^2 >= 0).
@raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}} or B{C{alphap2}} or no convergence.
@note: The arguments must satisfy I{B{C{k2}} + B{C{kp2}} = 1} and I{B{C{alpha2}} + B{C{alphap2}} = 1}. No checking is done that these conditions are met to enable accuracy to be maintained, e.g., when C{k} is very close to unity. ''' raise EllipticError('%s invalid: %r' % ('k2', k2))
raise EllipticError('%s invalid: %r' % ('alpha2', alpha2))
elif kp2 < 0: raise EllipticError('%s invalid: %r' % ('kp2', kp2))
elif alphap2 < 0: raise EllipticError('%s invalid: %r' % ('alphap2', alphap2))
# Values of complete elliptic integrals for k = 0,1 and alpha = 0,1 # K E D # k = 0: pi/2 pi/2 pi/4 # k = 1: inf 1 inf # Pi G H # k = 0, alpha = 0: pi/2 pi/2 pi/4 # k = 1, alpha = 0: inf 1 1 # k = 0, alpha = 1: inf inf pi/2 # k = 1, alpha = 1: inf inf inf # # G(0, k) = Pi(0, k) = H(1, k) = E(k) # H(0, k) = K(k) - D(k) # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2)) # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1)) # Pi(alpha2, 1) = inf # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2) # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3 # <https://DLMF.NIST.gov/19.25.E1> # Complete elliptic integral E(k), Carlson eq. 4.2 # <https://DLMF.NIST.gov/19.25.E1> # Complete elliptic integral K(k), Carlson eq. 4.1 # <https://DLMF.NIST.gov/19.25.E1> else: else:
# <https://DLMF.NIST.gov/19.25.E2> if alphap2: rj_3 = _RJ_3(0, kp2, 1, alphap2) # G(alpha2, k) self._cG = self.cK + (alpha2 - k2) * rj_3 # H(alpha2, k) self._cH = self.cK - alphap2 * rj_3 # Pi(alpha2, k) self._cPi = self.cK + alpha2 * rj_3 else: self._cG = self._cH = self._cPi = INF # XXX or NAN? self._cG = self._cH = _RC(1, alphap2) self._cPi = INF # XXX or NAN? else: # cH = cK - cD but this involves large cancellations # if k2 is close to 1. So write (for alpha2 = 0) # cH = int(cos(phi)**2/sqrt(1-k2*sin(phi)**2),phi,0,pi/2) # = 1/sqrt(1-k2) * int(sin(phi)**2/sqrt(1-k2/kp2*sin(phi)**2,...) # = 1/kp * D(i*k/kp) # and use D(k) = RD(0, kp2, 1) / 3 # so cH = 1/kp * RD(0, 1/kp2, 1) / 3 # = kp2 * RD(0, 1, kp2) / 3 # using <https://DLMF.NIST.gov/19.20.E18> # Equivalently # RF(x, 1) - RD(0, x, 1)/3 = x * RD(0, 1, x)/3 for x > 0 # For k2 = 1 and alpha2 = 0, we have # cH = int(cos(phi),...) = 1
'''The Jacobi elliptic function.
@param x: The argument (C{float}).
@return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with C{sn(B{x}, k), cn(B{x}, k), dn(B{x}, k)}.
@raise EllipticError: No convergence. ''' # Bulirsch's sncndn routine, p 89. if mc < 0: # PYCHOK no cover d = 1.0 - mc mc = -mc / d # /= -d chokes PyChecker d = sqrt(d) x *= d else: # This converges quadratically, max 6 trips else: raise EllipticError('no %s%r convergence' % ('sncndn', (x,))) if d: # PYCHOK no cover cn, dn = dn, cn sn = sn / d # /= d chokes PyChecker else:
'''(INTERNAL) Helper for C{.fEinv} and C{._fXf}. '''
'''3-Tuple C{(sn, cn, dn)} all C{float}. '''
'''(INTERNAL) Horner form for C{_RD} and C{_RJ} below. ''' # Polynomial is <https://DLMF.NIST.gov/19.36.E2> # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26 # - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20 + 45*E2**2*E3/272 # - 9*(E3*E4+E2*E5)/68) # converted to Horner form ...
'''(INTERNAL) Helper for C{_RD}, C{_RF} and C{_RJ}. '''
'''Degenerate symmetric integral of the first kind C{_RC}.
@return: C{_RC(x, y) = _RF(x, y, y)}.
@see: U{C{_RC} definition<https://DLMF.NIST.gov/19.2.E17>}. ''' # Defined only for y != 0 and x >= 0. # <https://DLMF.NIST.gov/19.2.E18> else: raise EllipticError('%s invalid: %s%r' % ('y', '_RC', (x, y)))
'''Degenerate symmetric integral of the third kind C{_RD} / 3. '''
'''Degenerate symmetric integral of the third kind C{_RD}.
@return: C{_RD(x, y, z) = _RJ(x, y, z, z)}.
@see: U{C{_RD} definition<https://DLMF.NIST.gov/19.16.E5>}. ''' # Carlson, eqs 2.28 - 2.34 else: raise EllipticError('no %s%r convergence' % ('_RD', (x, y, z)))
xy - 6 * z2, (xy * 3 - 8 * z2) * z, (xy - z2) * 3 * z2, xy * z2 * z)
'''Symmetric integral of the first kind C{_RF}.
@return: C{_RF(x, y)}.
@see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>}. ''' # Carlson, eqs 2.36 - 2.38
'''Symmetric integral of the first kind C{_RF}.
@return: C{_RF(x, y, z)}.
@see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>}. ''' # Carlson, eqs 2.2 - 2.7 else: raise EllipticError('no %s%r convergence' % ('_RF', (x, y, z)))
# Polynomial is <https://DLMF.NIST.gov/19.36.E1> # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16) # converted to Horner form ...
'''Symmetric integral of the second kind C{_RG}.
@return: C{_RG(x, y)}.
@see: U{C{_RG} definition<https://DLMF.NIST.gov/19.16.E3>} and in Carlson, eq. 1.5. ''' # Carlson, eqs 2.36 - 2.39
'''Symmetric integral of the second kind C{_RG}.
@return: C{_RG(x, y, z)}.
@see: U{C{_RG} definition<https://DLMF.NIST.gov/19.16.E3>} and in Carlson, eq. 1.5. ''' y, z = z, y # Carlson, eq 1.7 _RD_3(x, y, z) * (x - z) * (z - y), sqrt(x * y / z)) * 0.5
'''Symmetric integral of the third kind C{_RJ} / 3. '''
'''Symmetric integral of the third kind C{_RJ}.
@return: C{_RJ(x, y, z, p)}.
@see: U{C{_RJ} definition<https://DLMF.NIST.gov/19.16.E2>}. '''
# Carlson, eqs 2.17 - 2.25 else: raise EllipticError('no %s%r convergence' % ('_RJ', (x, y, z, p)))
fsum_(xyz, 2 * p * e2, 4 * p * p2), fsum_(xyz * 2, p * e2, 3 * p * p2) * p, p2 * xyz)
'''(INTERNAL) Helper for C{_RD}, C{_RF} and C{_RJ}. '''
# **) MIT License # # Copyright (C) 2016-2020 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |