Coverage for pygeodesy/elliptic.py : 97%

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 I{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>} (also available U{here <https://ArXiv.org/pdf/math/9409227.pdf>}), 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))
_no_convergence_, _0_0, _0_5, \ _1_0, _2_0, _3_0, _4_0, _5_0, _360_0 # from pygeodesy.streprs import unstr
sin, sqrt, tanh
'''Elliptic integral, function, convergence or other L{Elliptic} issue. '''
'''3-Tuple C{(sn, cn, dn)} all C{scalar}. '''
'''Elliptic integrals and functions.
@see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/ html/classGeographicLib_1_1EllipticFunction.html#details>}. '''
'''Constructor, specifying the C{modulus} and C{parameter}.
@kwarg k2: Modulus squared (C{float}, 0 <= k^2 <= 1). @kwarg alpha2: Parameter squared (C{float}, 0 <= α^2 <= 1). @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0). @kwarg 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). '''
'''Get α^2, the square of the parameter (C{float}). '''
'''Get α'^2, the square of the complementary parameter (C{float}). '''
'''Get Jahnke's complete integral C{D(k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E6>}. '''
'''Get the complete integral of the second kind C{E(k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}. '''
'''Get Legendre's complete geodesic longitude integral C{G(α^2, k)} (C{float}). '''
'''Get Cayley's complete geodesic longitude difference integral C{H(α^2, k)} (C{float}). '''
'''Get the complete integral of the first kind C{K(k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}. '''
'''Get the difference between the complete integrals of the first and second kinds, C{K(k) − E(k)} (C{float}). '''
'''Get the complete integral of the third kind C{Pi(α^2, k)} (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}. '''
'''The periodic Jahnke's incomplete elliptic integral.
@arg sn: sin(φ). @arg cn: cos(φ). @arg 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.
@arg sn: sin(φ). @arg cn: cos(φ). @arg 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.
@arg stau: sin(τ) @arg 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.
@arg sn: sin(φ). @arg cn: cos(φ). @arg 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.
@arg sn: sin(φ). @arg cn: cos(φ). @arg 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.
@arg sn: sin(φ). @arg cn: cos(φ). @arg 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.
@arg sn: sin(φ). @arg cn: cos(φ). @arg 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 _invokationError(t, sn, cn, dn)
cn, sn = -cn, -sn
'''Get epsilon (C{float}). '''
'''Jahnke's incomplete elliptic integral in terms of Jacobi elliptic functions.
@arg phi_or_sn: φ or sin(φ). @kwarg cn: C{None} or cos(φ). @kwarg 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.
@arg sn: sin(φ). @arg 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.
@arg phi_or_sn: φ or sin(φ). @kwarg cn: C{None} or cos(φ). @kwarg 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. ''' 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.
@arg 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.
@arg 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 _convergenceError(self.fEinv, x)
'''The incomplete integral of the first kind in terms of Jacobi elliptic functions.
@arg phi_or_sn: φ or sin(φ). @kwarg cn: C{None} or cos(φ). @kwarg 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.
@arg phi_or_sn: φ or sin(φ). @kwarg cn: C{None} or cos(φ). @kwarg 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.
@arg phi_or_sn: φ or sin(φ). @kwarg cn: C{None} or cos(φ). @kwarg 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.
@arg phi_or_sn: φ or sin(φ). @kwarg cn: C{None} or cos(φ). @kwarg 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 _invokationError(t, sn, cn, dn)
xi = 2 * cX - fX(sn, cn, dn) else:
'''Get the most recent C{Elliptic.fEinv} or C{Elliptic.sncndn} iteration number (C{int} or C{0} if not available/applicable). ''' return self._iteration
'''Get k^2, the square of the modulus (C{float}). '''
'''Get k'^2, the square of the complementary modulus (C{float}). '''
'''Reset the modulus and parameter.
@kwarg k2: modulus squared (C{float}, -INF <= k^2<= 1). @kwarg alpha2: parameter (C{float}, -INF <= α^2 <= 1). @kwarg kp2: complementary modulus squared (C{float}, k'^2 >= 0). @kwarg 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. '''
name='kp2', Error=EllipticError)
name='alphap2', Error=EllipticError)
# 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> # G(alpha2, k) # H(alpha2, k) # Pi(alpha2, k) else: self._cG = self._cH = _RC(1, alphap2) self._cPi = INF # XXX or NAN? else: 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.
@arg x: The argument (C{float}).
@return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with C{*n}C{(}B{C{x}}C{, k}C{)}.
@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 _convergenceError(self.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}. '''
def _convergenceError(where, *args): # PYCHOK no cover '''(INTERNAL) Return an L{EllipticError}. ''' return EllipticError(_no_convergence_, txt='%s%r' % (where.__name__, args)) # unstr
'''(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 ...
def _invokationError(name, *args): # PYCHOK no cover '''(INTERNAL) Return an L{EllipticError}. ''' return EllipticError('%s %s%r' % ('invokation', name, args)) # unstr
'''(INTERNAL) Helper for C{_RD}, C{_RF} and C{_RJ}. '''
'''Degenerate symmetric integral of the first kind C{_RC(x, y)}.
@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> elif y < 0: # <https://DLMF.NIST.gov/19.2.E20> r = asinh(sqrt(-x / y)) # atanh(sqrt(x / (x - y))) else: raise _invokationError(_RC.__name__, x, y)
'''Degenerate symmetric integral of the third kind C{_RD(x, y, z) / 3}. '''
'''Degenerate symmetric integral of the third kind C{_RD(x, y, z)}.
@return: C{_RD(x, y, z) = _RJ(x, y, z, z)}.
@see: U{C{_RD} definition<https://DLMF.NIST.gov/19.16.E5>} and U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. ''' # Carlson, eqs 2.28 - 2.34 else: raise _convergenceError(_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_(x, y)}.
@return: C{_RF(x, y)}.
@see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>}. ''' # Carlson, eqs 2.36 - 2.38
raise _convergenceError(_RF_, x, y)
'''Symmetric integral of the first kind C{_RF(x, y, z)}.
@return: C{_RF(x, y, z)}.
@see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>} and U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. ''' # Carlson, eqs 2.2 - 2.7 else: raise _convergenceError(_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_(x, y)}.
@return: C{_RG(x, y)}.
@see: U{C{_RG} definition<https://DLMF.NIST.gov/19.16.E3>} and U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. ''' # Carlson, eqs 2.36 - 2.39
raise _convergenceError(_RG_, x, y)
'''Symmetric integral of the second kind C{_RG(x, y, z)}.
@return: C{_RG(x, y, z)}.
@see: U{C{_RG} definition<https://DLMF.NIST.gov/19.16.E3>} and U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. ''' 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(x, y, z, p) / 3}. '''
'''Symmetric integral of the third kind C{_RJ(x, y, z, p)}.
@return: C{_RJ(x, y, z, p)}.
@see: U{C{_RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. '''
# Carlson, eqs 2.17 - 2.25 else: raise _convergenceError(_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. |