Coverage for pygeodesy/elliptic.py : 88%

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 the 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{alpha2}. The modulus is always passed as its square 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 as _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} k^2 <= 1). @keyword alpha2: Parameter squared (C{scalar} α^2 <= 1). @keyword kp2: Complementary modulus squared (C{float} k'^2 >= 0). @keyword alphap2: Complementary parameter α'2 (C{float} α'^2 >= 0).
@note: If only elliptic integrals of the first and second kinds are needed, then set B{C{alpha2}} = 0 (the default value); in this case, we have Π(φ, 0, k) = F(φ, k), G(φ, 0, k) = E(φ, k), and H(φ, 0, k) = F(φ, k) - D(φ, k).
@see: Method L{reset}. '''
def alpha2(self): '''Get the parameter B{C{alpha2}} (C{float}). '''
def alphap2(self): '''Get the complementary parameter B{C{alphap2}} (C{float}). '''
def cD(self): '''Get Jahnke's complete integral C{D(k)}, U{defined<https://DLMF.NIST.gov/19.2.E6>}. ''' return self._Dc
def cE(self): '''Get the complete integral of the second kind C{E(k)}, U{defined<https://DLMF.NIST.gov/19.2.E5>}. '''
def cG(self): '''Get Legendre's complete geodesic longitude integral C{G(alpha2, k)}. '''
def cH(self): '''Get Cayley's complete geodesic longitude difference integral C{H(alpha2, k)}. '''
def cK(self): '''Get the complete integral of the first kind C{K(k)}, 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)}. '''
def cPi(self): '''Get the complete integral of the third kind C{Pi(alpha2, k)}, 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 sin2φ).
@return: Periodic function π D(φ, k) / (2 D(k)) - φ. ''' return self._deltaX(sn, cn, dn, self.cD, self.fD)
'''The periodic incomplete integral of the second kind.
@param sn: sinφ. @param cn: cosφ. @param dn: sqrt(1 − k2 sin2φ).
@return: Periodic function 𝜋π E(φ, k) / (2 E(k)) - φ. ''' return self._deltaX(sn, cn, dn, self.cE, self.fE)
'''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) - τ. ''' # Function is periodic with period pi t = atan2(-stau, -ctau) if ctau < 0 else atan2(stau, ctau) return self.fEinv(t * self._Ec / PI_2) - t
'''The periodic incomplete integral of the first kind.
@param sn: sinφ. @param cn: cosφ. @param dn: sqrt(1 − k2 sin2φ).
@return: Periodic function π F(φ, k) / (2 K(k)) - φ. ''' return self._deltaX(sn, cn, dn, self.cK, self.fF)
'''Legendre's periodic geodesic longitude integral.
@param sn: sinφ. @param cn: cosφ. @param dn: sqrt(1 − k2 sin2φ).
@return: Periodic function π G(φ, k) / (2 G(k)) - φ. ''' return self._deltaX(sn, cn, dn, self.cG, self.fG)
'''Cayley's periodic geodesic longitude difference integral.
@param sn: sinφ. @param cn: cosφ. @param dn: sqrt(1 − k2 sin2φ).
@return: Periodic function π H(φ, k) / (2 H(k)) - φ. ''' return self._deltaX(sn, cn, dn, self.cH, self.fH)
'''The periodic incomplete integral of the third kind.
@param sn: sinφ. @param cn: cosφ. @param dn: sqrt(1 − k2 sin2φ).
@return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ. ''' return self._deltaX(sn, cn, dn, self.cPi, self.fPi)
'''(INTERNAL) Helper for C{.deltaD} thru C{.deltaPi}. ''' if cn < 0: cn, sn = -cn, -sn return fX(sn, cn, dn) * PI_2 / cX - atan2(sn, cn)
'''Jahnke's incomplete elliptic integral in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sinφ. @keyword cn: cosφ. @keyword dn: sqrt(1 − k2 sin2φ).
@return: D(φ, k) as though φ ∈ (−π, π], U{defined<https://DLMF.NIST.gov/19.2.E4>}. ''' def _fD(sn, cn, dn): return abs(sn) * sn**2 * _RD_3(cn**2, dn**2, 1)
return self._fXf(phi_or_sn, cn, dn, self.cD, self.deltaD, _fD)
'''The C{Delta} amplitude function.
@param sn: sinφ. @param cn: cosφ.
@return: C{sqrt(1 − k2 sin2φ)}. ''' (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: cosφ. @keyword dn: sqrt(1 − k2 sin2φ).
@return: E(φ, k) as though φ ∈ (−π, π], U{defined<https://DLMF.NIST.gov/19.2.E5>}. ''' # 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(π deg/180, k). '''
'''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}}. ''' # 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 convergence' % ('fEinv',))
'''The incomplete integral of the first kind in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sinφ. @keyword cn: cosφ. @keyword dn: sqrt(1 − k2 sin2φ).
@return: F(φ, k) as though φ ∈ (−π, π], U{defined<https://DLMF.NIST.gov/19.2.E4>}. '''
self.deltaF, _fF)
'''Legendre's geodesic longitude integral in terms of Jacobi elliptic functions.
@param phi_or_sn: φ or sinφ. @keyword cn: cosφ. @keyword dn: sqrt(1 − k2 sin2φ).
@return: G(φ, k) as though φ ∈ (−π, π].
@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: cosφ. @keyword dn: sqrt(1 − k2 sin2φ).
@return: H(φ, k) as though φ ∈ (−π, π].
@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: cosφ. @keyword dn: sqrt(1 − k2 sin2φ).
@return: Π(φ, α2, k) as though φ ∈ (−π, π]. ''' 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 n = self.classname + '.f' + deltaX.__name__[5:] raise EllipticError('%s invalid: %s%r' % ('args', n, (sn, cn, dn)))
elif cn < 0: xi = 2 * cX - fX(sn, cn, dn) else: xi = cX
def k2(self): '''Get the square of the modulus (C{float}). '''
def kp2(self): '''Get the square of the complementary modulus (C{float}). ''' return self._kp2
'''Reset the modulus and parameter.
@keyword k2: modulus squared (C{float} k2 <= 1). @keyword alpha2: parameter (C{float} α2 <= 1). @keyword kp2: complementary modulus squared (C{float} k'2 >= 0). @keyword alphap2: complementary parameter α'2 (C{float} α'2 >= 0).
@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)) else:
raise EllipticError('%s invalid: %r' % ('alpha2', alpha2)) else:
elif kp2 < 0: raise EllipticError('%s invalid: %r' % ('kp2', kp2)) else: self._kp2 = kp2 = float(kp2)
elif alphap2 < 0: raise EllipticError('%s invalid: %r' % ('alphap2', alphap2)) else: self._alphap2 = alphap2 = float(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) = E(k) # H( 0, k) = K(k) - D(k) # Pi(0, k) = K(k) # H( 1, k) = K(k) # Pi(alpha2, 0) = pi/(2*sqrt(1-alpha2)) # G( alpha2, 0) = pi/(2*sqrt(1-alpha2)) # H( alpha2, 0) = pi/(2*(1 + sqrt(1-alpha2))) # Pi(alpha2, 1) = inf # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2) # D(k) = (K(k) - E(k))/k^2, 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(alpha^2, k) self._Gc = self._Kc + (alpha2 - k2) * rj_3 # H(alpha^2, k) self._Hc = self._Kc - alphap2 * rj_3 # Pi(alpha^2, k) self._Pic = self._Kc + alpha2 * rj_3 else: self._Gc = self._Hc = self._Pic = _INF # XXX _NAN? self._Gc = self._Hc = _RC(1, alphap2) self._Pic = _INF # XXX _NAN? else: # Hc = Kc - Dc but this involves large cancellations # if k2 is close to 1. So write (for alpha2 = 0) # Hc = 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 Hc = 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 # Hc = int(cos(phi),...) = 1
'''The Jacobi elliptic function.
@param x: The argument (C{float}).
@return: 3-Tuple C{(sn(x, k), cn(x, k), dn(x, k))}. ''' # Bulirsch's sncndn routine, p 89. if mc < 0: # PYCHOK no cover d = 1.0 - mc mc /= -d d = sqrt(d) x *= d else: # This converges quadratically, max 6 trips else: raise EllipticError('no %s convergence' % ('sncndn',)) cn, dn = dn, cn sn /= d else:
'''(INTERNAL) Helper for C{.fEinv} and C{._fXf}. '''
'''(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) fsum_(612612 * e2, -540540 * e3, -556920) * e4, fsum_(306306 * e3, 675675 * e2**2, -706860 * e2, 680680) * e3, fsum_(417690 * e2, -255255 * e2**2, -875160) * e2, 4084080) / (4084080 * e1) + e0
'''(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: %r' % ('y', (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 convergence' % ('RD',))
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, z)}.
@see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>}. ''' # Carlson, eqs 2.36 - 2.38
# Carlson, eqs 2.2 - 2.7 else: raise EllipticError('no %s convergence' % ('RF',))
# 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) # convert to Horner form ... fsum_(10010 * e2, -5775 * e2**2, -24024) * e2, 240240) / (240240 * sqrt(T[0]))
'''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. ''' # Carlson, eqs 2.36 - 2.39
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 convergence' % ('RJ',))
e2, 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. |