Coverage for pygeodesy/elliptic.py: 99%
437 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
2# -*- coding: utf-8 -*-
4u'''I{Karney}'s elliptic functions and integrals.
6Class L{Elliptic} transcoded from I{Charles Karney}'s C++ class U{EllipticFunction
7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1EllipticFunction.html>}
8to pure Python, including symmetric integrals L{Elliptic.fRC}, L{Elliptic.fRD},
9L{Elliptic.fRF}, L{Elliptic.fRG} and L{Elliptic.fRJ} as C{static methods}.
11Python method names follow the C++ member functions, I{except}:
13 - member functions I{without arguments} are mapped to Python properties
14 prefixed with C{"c"}, for example C{E()} is property C{cE},
16 - member functions with 1 or 3 arguments are renamed to Python methods
17 starting with an C{"f"}, example C{E(psi)} to C{fE(psi)} and C{E(sn,
18 cn, dn)} to C{fE(sn, cn, dn)},
20 - other Python method names conventionally start with a lower-case
21 letter or an underscore if private.
23Following is a copy of I{Karney}'s U{EllipticFunction.hpp
24<https://GeographicLib.SourceForge.io/C++/doc/EllipticFunction_8hpp_source.html>}
25file C{Header}.
27Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023)
28and licensed under the MIT/X11 License. For more information, see the
29U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
31B{Elliptic integrals and functions.}
33This provides the elliptic functions and integrals needed for
34C{Ellipsoid}, C{GeodesicExact}, and C{TransverseMercatorExact}. Two
35categories of function are provided:
37 - functions to compute U{symmetric elliptic integrals
38 <https://DLMF.NIST.gov/19.16.i>}
40 - methods to compute U{Legrendre's elliptic integrals
41 <https://DLMF.NIST.gov/19.2.ii>} and U{Jacobi elliptic
42 functions<https://DLMF.NIST.gov/22.2>}.
44In the latter case, an object is constructed giving the modulus
45C{k} (and optionally the parameter C{alpha}). The modulus (and
46parameter) are always passed as squares which allows C{k} to be
47pure imaginary. (Confusingly, Abramowitz and Stegun call C{m = k**2}
48the "parameter" and C{n = alpha**2} the "characteristic".)
50In geodesic applications, it is convenient to separate the incomplete
51integrals into secular and periodic components, e.g.
53I{C{E(phi, k) = (2 E(k) / pi) [ phi + delta E(phi, k) ]}}
55where I{C{delta E(phi, k)}} is an odd periodic function with
56period I{C{pi}}.
58The computation of the elliptic integrals uses the algorithms given
59in U{B. C. Carlson, Computation of real or complex elliptic integrals
60<https://DOI.org/10.1007/BF02198293>} (also available U{here
61<https://ArXiv.org/pdf/math/9409227.pdf>}), Numerical Algorithms 10,
6213--26 (1995) with the additional optimizations given U{here
63<https://DLMF.NIST.gov/19.36.i>}.
65The computation of the Jacobi elliptic functions uses the algorithm
66given in U{R. Bulirsch, Numerical Calculation of Elliptic Integrals
67and Elliptic Functions<https://DOI.org/10.1007/BF01397975>},
68Numerische Mathematik 7, 78--90 (1965).
70The notation follows U{NIST Digital Library of Mathematical Functions
71<https://DLMF.NIST.gov>} chapters U{19<https://DLMF.NIST.gov/19>} and
72U{22<https://DLMF.NIST.gov/22>}.
73'''
74# make sure int/int division yields float quotient, see .basics
75from __future__ import division as _; del _ # PYCHOK semicolon
77from pygeodesy.basics import copysign0, map2, neg
78from pygeodesy.constants import EPS, INF, PI, PI_2, PI_4, \
79 _EPStol as _TolJAC, _over, \
80 _0_0, _0_125, _0_25, _0_5, _1_64th, \
81 _1_0, _2_0, _N_2_0, _3_0, _4_0, _6_0, \
82 _8_0, _180_0, _360_0
83from pygeodesy.errors import _ValueError, _xkwds_pop
84from pygeodesy.fmath import fdot, hypot1
85from pygeodesy.fsums import Fsum, _sum
86from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _SPACE_
87from pygeodesy.karney import _ALL_LAZY, _signBit
88# from pygeodesy.lazily import _ALL_LAZY # from .karney
89from pygeodesy.named import _Named, _NamedTuple
90from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all
91from pygeodesy.streprs import Fmt, unstr
92from pygeodesy.units import Scalar, Scalar_
93from pygeodesy.utily import sincos2, sincos2d
95from math import asinh, atan, atan2, ceil, cosh, fabs, floor, sin, sqrt, tanh
97__all__ = _ALL_LAZY.elliptic
98__version__ = '23.08.28'
100_invokation_ = 'invokation'
101_TolRD = pow(EPS * 0.002, _0_125) # 8th root: quadquadratic?, qqrt, oqrt?
102_TolRF = pow(EPS * 0.030, _0_125) # 4th root: biquadratic, bqrt?
103_TolRG0 = _TolJAC * 2.7
104_TRIPS = 21 # Max depth, 7 might be sufficient
107class _CIs(object):
108 '''(INTERAL) Hold the complete integrals.
109 '''
110 def __init__(self, **kwds):
111 self.__dict__ = kwds
114class _Dsum(list):
115 '''(INTERNAL) Deferred C{Fsum}.
116 '''
117 def __call__(self, s):
118 return Fsum(*self).fmul(s)
120 def __iadd__(self, x):
121 list.append(self, x)
122 return self
125class Elliptic(_Named):
126 '''Elliptic integrals and functions.
128 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/
129 html/classGeographicLib_1_1EllipticFunction.html#details>}.
130 '''
131# _alpha2 = 0
132# _alphap2 = 0
133# _eps = EPS
134# _k2 = 0
135# _kp2 = 0
137 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None, name=NN):
138 '''Constructor, specifying the C{modulus} and C{parameter}.
140 @kwarg name: Optional name (C{str}).
142 @see: Method L{Elliptic.reset} for further details.
144 @note: If only elliptic integrals of the first and second kinds
145 are needed, use C{B{alpha2}=0}, the default value. In
146 that case, we have C{Π(φ, 0, k) = F(φ, k), G(φ, 0, k) =
147 E(φ, k)} and C{H(φ, 0, k) = F(φ, k) - D(φ, k)}.
148 '''
149 self.reset(k2=k2, alpha2=alpha2, kp2=kp2, alphap2=alphap2)
151 if name:
152 self.name = name
154 @Property_RO
155 def alpha2(self):
156 '''Get α^2, the square of the parameter (C{float}).
157 '''
158 return self._alpha2
160 @Property_RO
161 def alphap2(self):
162 '''Get α'^2, the square of the complementary parameter (C{float}).
163 '''
164 return self._alphap2
166 @Property_RO
167 def cD(self):
168 '''Get Jahnke's complete integral C{D(k)} (C{float}),
169 U{defined<https://DLMF.NIST.gov/19.2.E6>}.
170 '''
171 return self._reset_cDcEcKcKEeps.cD
173 @Property_RO
174 def cE(self):
175 '''Get the complete integral of the second kind C{E(k)}
176 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}.
177 '''
178 return self._reset_cDcEcKcKEeps.cE
180 @Property_RO
181 def cG(self):
182 '''Get Legendre's complete geodesic longitude integral
183 C{G(α^2, k)} (C{float}).
184 '''
185 return self._reset_cGcHcPi.cG
187 @Property_RO
188 def cH(self):
189 '''Get Cayley's complete geodesic longitude difference integral
190 C{H(α^2, k)} (C{float}).
191 '''
192 return self._reset_cGcHcPi.cH
194 @Property_RO
195 def cK(self):
196 '''Get the complete integral of the first kind C{K(k)}
197 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}.
198 '''
199 return self._reset_cDcEcKcKEeps.cK
201 @Property_RO
202 def cKE(self):
203 '''Get the difference between the complete integrals of the
204 first and second kinds, C{K(k) − E(k)} (C{float}).
205 '''
206 return self._reset_cDcEcKcKEeps.cKE
208 @Property_RO
209 def cPi(self):
210 '''Get the complete integral of the third kind C{Pi(α^2, k)}
211 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}.
212 '''
213 return self._reset_cGcHcPi.cPi
215 def deltaD(self, sn, cn, dn):
216 '''The periodic Jahnke's incomplete elliptic integral.
218 @arg sn: sin(φ).
219 @arg cn: cos(φ).
220 @arg dn: sqrt(1 − k2 * sin(2φ)).
222 @return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}).
224 @raise EllipticError: Invalid invokation or no convergence.
225 '''
226 return _deltaX(sn, cn, dn, self.cD, self.fD)
228 def deltaE(self, sn, cn, dn):
229 '''The periodic incomplete integral of the second kind.
231 @arg sn: sin(φ).
232 @arg cn: cos(φ).
233 @arg dn: sqrt(1 − k2 * sin(2φ)).
235 @return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}).
237 @raise EllipticError: Invalid invokation or no convergence.
238 '''
239 return _deltaX(sn, cn, dn, self.cE, self.fE)
241 def deltaEinv(self, stau, ctau):
242 '''The periodic inverse of the incomplete integral of the second kind.
244 @arg stau: sin(τ)
245 @arg ctau: cos(τ)
247 @return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}).
249 @raise EllipticError: No convergence.
250 '''
251 # Function is periodic with period pi
252 t = atan2(-stau, -ctau) if _signBit(ctau) else atan2(stau, ctau)
253 return self.fEinv(t * self.cE / PI_2) - t
255 def deltaF(self, sn, cn, dn):
256 '''The periodic incomplete integral of the first kind.
258 @arg sn: sin(φ).
259 @arg cn: cos(φ).
260 @arg dn: sqrt(1 − k2 * sin(2φ)).
262 @return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}).
264 @raise EllipticError: Invalid invokation or no convergence.
265 '''
266 return _deltaX(sn, cn, dn, self.cK, self.fF)
268 def deltaG(self, sn, cn, dn):
269 '''Legendre's periodic geodesic longitude integral.
271 @arg sn: sin(φ).
272 @arg cn: cos(φ).
273 @arg dn: sqrt(1 − k2 * sin(2φ)).
275 @return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}).
277 @raise EllipticError: Invalid invokation or no convergence.
278 '''
279 return _deltaX(sn, cn, dn, self.cG, self.fG)
281 def deltaH(self, sn, cn, dn):
282 '''Cayley's periodic geodesic longitude difference integral.
284 @arg sn: sin(φ).
285 @arg cn: cos(φ).
286 @arg dn: sqrt(1 − k2 * sin(2φ)).
288 @return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}).
290 @raise EllipticError: Invalid invokation or no convergence.
291 '''
292 return _deltaX(sn, cn, dn, self.cH, self.fH)
294 def deltaPi(self, sn, cn, dn):
295 '''The periodic incomplete integral of the third kind.
297 @arg sn: sin(φ).
298 @arg cn: cos(φ).
299 @arg dn: sqrt(1 − k2 * sin(2φ)).
301 @return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ
302 (C{float}).
304 @raise EllipticError: Invalid invokation or no convergence.
305 '''
306 return _deltaX(sn, cn, dn, self.cPi, self.fPi)
308 @Property_RO
309 def eps(self):
310 '''Get epsilon (C{float}).
311 '''
312 return self._reset_cDcEcKcKEeps.eps
314 def fD(self, phi_or_sn, cn=None, dn=None):
315 '''Jahnke's incomplete elliptic integral in terms of
316 Jacobi elliptic functions.
318 @arg phi_or_sn: φ or sin(φ).
319 @kwarg cn: C{None} or cos(φ).
320 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
322 @return: D(φ, k) as though φ ∈ (−π, π] (C{float}),
323 U{defined<https://DLMF.NIST.gov/19.2.E4>}.
325 @raise EllipticError: Invalid invokation or no convergence.
326 '''
327 def _fD(sn, cn, dn):
328 r = fabs(sn)**3
329 if r:
330 r = float(_RD(self, cn**2, dn**2, _1_0, _3_0 / r))
331 return r
333 return self._fXf(phi_or_sn, cn, dn, self.cD,
334 self.deltaD, _fD)
336 def fDelta(self, sn, cn):
337 '''The C{Delta} amplitude function.
339 @arg sn: sin(φ).
340 @arg cn: cos(φ).
342 @return: sqrt(1 − k2 * sin(2φ)) (C{float}).
343 '''
344 k2 = self.k2
345 s = (_1_0 - sn**2 * k2) if k2 < 0 else (self.kp2
346 + ((cn**2 * k2) if k2 > 0 else _0_0))
347 return sqrt(s) if s else _0_0
349 def fE(self, phi_or_sn, cn=None, dn=None):
350 '''The incomplete integral of the second kind in terms of
351 Jacobi elliptic functions.
353 @arg phi_or_sn: φ or sin(φ).
354 @kwarg cn: C{None} or cos(φ).
355 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
357 @return: E(φ, k) as though φ ∈ (−π, π] (C{float}),
358 U{defined<https://DLMF.NIST.gov/19.2.E5>}.
360 @raise EllipticError: Invalid invokation or no convergence.
361 '''
362 def _fE(sn, cn, dn):
363 '''(INTERNAL) Core of C{.fE}.
364 '''
365 if sn:
366 sn2, cn2, dn2 = sn**2, cn**2, dn**2
367 kp2, k2 = self.kp2, self.k2
368 if k2 <= 0: # Carlson, eq. 4.6, <https://DLMF.NIST.gov/19.25.E9>
369 Ei = _RF3(self, cn2, dn2, _1_0)
370 if k2:
371 Ei -= _RD(self, cn2, dn2, _1_0, _3over(k2, sn2))
372 elif kp2 >= 0: # k2 > 0, <https://DLMF.NIST.gov/19.25.E10>
373 Ei = _over(k2 * fabs(cn), dn) # float
374 if kp2:
375 Ei += (_RD( self, cn2, _1_0, dn2, _3over(k2, sn2)) +
376 _RF3(self, cn2, dn2, _1_0)) * kp2
377 else: # kp2 < 0, <https://DLMF.NIST.gov/19.25.E11>
378 Ei = _over(dn, fabs(cn))
379 Ei -= _RD(self, dn2, _1_0, cn2, _3over(kp2, sn2))
380 Ei *= fabs(sn)
381 ei = float(Ei)
382 else: # PYCHOK no cover
383 ei = _0_0
384 return ei
386 return self._fXf(phi_or_sn, cn, dn, self.cE,
387 self.deltaE, _fE)
389 def fEd(self, deg):
390 '''The incomplete integral of the second kind with
391 the argument given in degrees.
393 @arg deg: Angle (C{degrees}).
395 @return: E(π B{C{deg}}/180, k) (C{float}).
397 @raise EllipticError: No convergence.
398 '''
399 if fabs(deg) < _180_0:
400 e = _0_0
401 else: # PYCHOK no cover
402 e = ceil(deg / _360_0 - _0_5)
403 deg -= e * _360_0
404 e *= self.cE * _4_0
405 sn, cn = sincos2d(deg)
406 return self.fE(sn, cn, self.fDelta(sn, cn)) + e
408 def fEinv(self, x):
409 '''The inverse of the incomplete integral of the second kind.
411 @arg x: Argument (C{float}).
413 @return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}}
414 (C{float}).
416 @raise EllipticError: No convergence.
417 '''
418 E2 = self.cE * _2_0
419 n = floor(x / E2 + _0_5)
420 r = x - E2 * n # r in [-cE, cE)
421 # linear approximation
422 phi = PI * r / E2 # phi in [-PI_2, PI_2)
423 Phi = Fsum(phi)
424 # first order correction
425 phi = Phi.fsum_(self.eps * sin(phi * _2_0) / _N_2_0)
426 # For kp2 close to zero use asin(r / cE) or J. P. Boyd,
427 # Applied Math. and Computation 218, 7005-7013 (2012)
428 # <https://DOI.org/10.1016/j.amc.2011.12.021>
429 _Phi2, self._iteration = Phi.fsum2_, 0 # aggregate
430 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC
431 sn, cn, dn = self._sncndnPhi(phi)
432 if sn:
433 sn = self.fE(sn, cn, dn)
434 phi, e = _Phi2((r - sn) / dn)
435 if fabs(e) < _TolJAC:
436 _iterations(self, i)
437 break
438 else: # PYCHOK no cover
439 raise _convergenceError(e, _TolJAC, self.fEinv, x)
440 return Phi.fsum_(n * PI) if n else phi
442 def fF(self, phi_or_sn, cn=None, dn=None):
443 '''The incomplete integral of the first kind in terms of
444 Jacobi elliptic functions.
446 @arg phi_or_sn: φ or sin(φ).
447 @kwarg cn: C{None} or cos(φ).
448 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
450 @return: F(φ, k) as though φ ∈ (−π, π] (C{float}),
451 U{defined<https://DLMF.NIST.gov/19.2.E4>}.
453 @raise EllipticError: Invalid invokation or no convergence.
454 '''
455 def _fF(sn, cn, dn):
456 r = fabs(sn)
457 if r:
458 r = float(_RF3(self, cn**2, dn**2, _1_0).fmul(r))
459 return r
461 return self._fXf(phi_or_sn, cn, dn, self.cK,
462 self.deltaF, _fF)
464 def fG(self, phi_or_sn, cn=None, dn=None):
465 '''Legendre's geodesic longitude integral in terms of
466 Jacobi elliptic functions.
468 @arg phi_or_sn: φ or sin(φ).
469 @kwarg cn: C{None} or cos(φ).
470 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
472 @return: G(φ, k) as though φ ∈ (−π, π] (C{float}).
474 @raise EllipticError: Invalid invokation or no convergence.
476 @note: Legendre expresses the longitude of a point on the
477 geodesic in terms of this combination of elliptic
478 integrals in U{Exercices de Calcul Intégral, Vol 1
479 (1811), p 181<https://Books.Google.com/books?id=
480 riIOAAAAQAAJ&pg=PA181>}.
482 @see: U{Geodesics in terms of elliptic integrals<https://
483 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>}
484 for the expression for the longitude in terms of this function.
485 '''
486 return self._fXa(phi_or_sn, cn, dn, self.alpha2 - self.k2,
487 self.cG, self.deltaG)
489 def fH(self, phi_or_sn, cn=None, dn=None):
490 '''Cayley's geodesic longitude difference integral in terms of
491 Jacobi elliptic functions.
493 @arg phi_or_sn: φ or sin(φ).
494 @kwarg cn: C{None} or cos(φ).
495 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
497 @return: H(φ, k) as though φ ∈ (−π, π] (C{float}).
499 @raise EllipticError: Invalid invokation or no convergence.
501 @note: Cayley expresses the longitude difference of a point
502 on the geodesic in terms of this combination of
503 elliptic integrals in U{Phil. Mag. B{40} (1870), p 333
504 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}.
506 @see: U{Geodesics in terms of elliptic integrals<https://
507 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>}
508 for the expression for the longitude in terms of this function.
509 '''
510 return self._fXa(phi_or_sn, cn, dn, -self.alphap2,
511 self.cH, self.deltaH)
513 def fPi(self, phi_or_sn, cn=None, dn=None):
514 '''The incomplete integral of the third kind in terms of
515 Jacobi elliptic functions.
517 @arg phi_or_sn: φ or sin(φ).
518 @kwarg cn: C{None} or cos(φ).
519 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
521 @return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}).
523 @raise EllipticError: Invalid invokation or no convergence.
524 '''
525 if dn is None and cn is not None: # and isscalar(phi_or_sn)
526 dn = self.fDelta(phi_or_sn, cn) # in .triaxial
527 return self._fXa(phi_or_sn, cn, dn, self.alpha2,
528 self.cPi, self.deltaPi)
530 def _fXa(self, phi_or_sn, cn, dn, aX, cX, _deltaX):
531 '''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}.
532 '''
533 def _fX(sn, cn, dn):
534 if sn:
535 cn2, dn2 = cn**2, dn**2
536 R = _RF3(self, cn2, dn2, _1_0)
537 if aX:
538 sn2 = sn**2
539 p = sn2 * self.alphap2 + cn2
540 R += _RJ(self, cn2, dn2, _1_0, p, _3over(aX, sn2))
541 r = float(R.fmul(fabs(sn)))
542 else: # PYCHOK no cover
543 r = _0_0
544 return r
546 return self._fXf(phi_or_sn, cn, dn, cX, _deltaX, _fX)
548 def _fXf(self, phi_or_sn, cn, dn, cX, _deltaX, _fX):
549 '''(INTERNAL) Helper for C{f.D}, C{.fE}, C{.fF} and C{._fXa}.
550 '''
551 self._iteration = 0 # aggregate
552 phi = sn = phi_or_sn
553 if cn is dn is None: # fX(phi) call
554 sn, cn, dn = self._sncndnPhi(phi)
555 if fabs(phi) >= PI:
556 if cX:
557 cX *= (_deltaX(sn, cn, dn) + phi) / PI_2
558 return cX
559 # fall through
560 elif cn is None or dn is None:
561 n = NN(_f_, _deltaX.__name__[5:])
562 raise _invokationError(n, sn, cn, dn)
564 if _signBit(cn): # enforce usual trig-like symmetries
565 xi = cX * _2_0 - _fX(sn, cn, dn)
566 elif cn > 0:
567 xi = _fX(sn, cn, dn)
568 else:
569 xi = cX
570 return copysign0(xi, sn)
572 @Property_RO
573 def k2(self):
574 '''Get k^2, the square of the modulus (C{float}).
575 '''
576 return self._k2
578 @Property_RO
579 def kp2(self):
580 '''Get k'^2, the square of the complementary modulus (C{float}).
581 '''
582 return self._kp2
584 def reset(self, k2=0, alpha2=0, kp2=None, alphap2=None): # MCCABE 13
585 '''Reset the modulus, parameter and the complementaries.
587 @kwarg k2: Modulus squared (C{float}, NINF <= k^2 <= 1).
588 @kwarg alpha2: Parameter squared (C{float}, NINF <= α^2 <= 1).
589 @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0).
590 @kwarg alphap2: Complementary parameter squared (C{float}, α'^2 >= 0).
592 @raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}}
593 or B{C{alphap2}}.
595 @note: The arguments must satisfy C{B{k2} + B{kp2} = 1} and
596 C{B{alpha2} + B{alphap2} = 1}. No checking is done
597 that these conditions are met to enable accuracy to be
598 maintained, e.g., when C{k} is very close to unity.
599 '''
600 if self.__dict__:
601 _update_all(self, _Named.iteration._uname, Base=Property_RO)
603 self._k2 = Scalar_(k2=k2, Error=EllipticError, low=None, high=_1_0)
604 self._kp2 = Scalar_(kp2=((_1_0 - k2) if kp2 is None else kp2), Error=EllipticError)
606 self._alpha2 = Scalar_(alpha2=alpha2, Error=EllipticError, low=None, high=_1_0)
607 self._alphap2 = Scalar_(alphap2=((_1_0 - alpha2) if alphap2 is None else alphap2),
608 Error=EllipticError)
610 # Values of complete elliptic integrals for k = 0,1 and alpha = 0,1
611 # K E D
612 # k = 0: pi/2 pi/2 pi/4
613 # k = 1: inf 1 inf
614 # Pi G H
615 # k = 0, alpha = 0: pi/2 pi/2 pi/4
616 # k = 1, alpha = 0: inf 1 1
617 # k = 0, alpha = 1: inf inf pi/2
618 # k = 1, alpha = 1: inf inf inf
619 #
620 # G(0, k) = Pi(0, k) = H(1, k) = E(k)
621 # H(0, k) = K(k) - D(k)
622 # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2))
623 # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1))
624 # Pi(alpha2, 1) = inf
625 # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2)
627 @Property_RO
628 def _reset_cDcEcKcKEeps(self):
629 '''(INTERNAL) Get the complete integrals D, E, K and KE plus C{eps}.
630 '''
631 k2 = self.k2
632 if k2:
633 kp2 = self.kp2
634 if kp2:
635 self._iteration = 0
636 # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3
637 # <https://DLMF.NIST.gov/19.25.E1>
638 D = _RD(self, _0_0, kp2, _1_0, _3_0)
639 cD = float(D)
640 # Complete elliptic integral E(k), Carlson eq. 4.2
641 # <https://DLMF.NIST.gov/19.25.E1>
642 cE = float(_RG2(self, kp2, _1_0))
643 # Complete elliptic integral K(k), Carlson eq. 4.1
644 # <https://DLMF.NIST.gov/19.25.E1>
645 cK = _rF2(self, kp2, _1_0)
646 cKE = float(D.fmul(k2))
647 eps = k2 / (sqrt(kp2) + _1_0)**2
648 else: # PYCHOK no cover
649 cD = cK = cKE = INF
650 cE = _1_0
651 eps = k2
652 else: # PYCHOK no cover
653 cD = PI_4
654 cE = cK = PI_2
655 cKE = _0_0 # k2 * cD
656 eps = EPS
658 return _CIs(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps)
660 @Property_RO
661 def _reset_cGcHcPi(self):
662 '''(INTERNAL) Get the complete integrals G, H and Pi.
663 '''
664 self._iteration = 0
665 alpha2 = self.alpha2
666 if alpha2:
667 alphap2 = self.alphap2
668 if alphap2:
669 kp2 = self.kp2
670 if kp2: # <https://DLMF.NIST.gov/19.25.E2>
671 cK = self.cK
672 Rj = _RJ(self, _0_0, kp2, _1_0, alphap2, _3_0)
673 cG = float(Rj * (alpha2 - self.k2) + cK) # G(alpha2, k)
674 cH = -float(Rj * alphap2 - cK) # H(alpha2, k)
675 cPi = float(Rj * alpha2 + cK) # Pi(alpha2, k)
676 else: # PYCHOK no cover
677 cG = cH = _rC(self, _1_0, alphap2)
678 cPi = INF # XXX or NAN?
679 else: # PYCHOK no cover
680 cG = cH = cPi = INF # XXX or NAN?
681 else:
682 cG, cPi, kp2 = self.cE, self.cK, self.kp2
683 # H = K - D but this involves large cancellations if k2 is near 1.
684 # So write (for alpha2 = 0)
685 # H = int(cos(phi)**2/sqrt(1-k2*sin(phi)**2),phi,0,pi/2)
686 # = 1/sqrt(1-k2) * int(sin(phi)**2/sqrt(1-k2/kp2*sin(phi)**2,...)
687 # = 1/kp * D(i*k/kp)
688 # and use D(k) = RD(0, kp2, 1) / 3
689 # so H = 1/kp * RD(0, 1/kp2, 1) / 3
690 # = kp2 * RD(0, 1, kp2) / 3
691 # using <https://DLMF.NIST.gov/19.20.E18>. Equivalently
692 # RF(x, 1) - RD(0, x, 1)/3 = x * RD(0, 1, x)/3 for x > 0
693 # For k2 = 1 and alpha2 = 0, we have
694 # H = int(cos(phi),...) = 1
695 cH = float(_RD(self, _0_0, _1_0, kp2, _3_0 / kp2)) if kp2 else _1_0
697 return _CIs(cG=cG, cH=cH, cPi=cPi)
699 def sncndn(self, x):
700 '''The Jacobi elliptic function.
702 @arg x: The argument (C{float}).
704 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with
705 C{*n(B{x}, k)}.
707 @raise EllipticError: No convergence.
708 '''
709 self._iteration = 0 # reset
710 # Bulirsch's sncndn routine, p 89.
711 if self.kp2:
712 c, d, cd, mn_ = self._sncndnBulirsch4
713 dn = _1_0
714 sn, cn = sincos2(x * cd)
715 if sn:
716 a = cn / sn
717 c *= a
718 for m, n in mn_:
719 a *= c
720 c *= dn
721 dn = (n + a) / (m + a)
722 a = c / m
723 sn = copysign0(_1_0 / hypot1(c), sn) # _signBit(sn)
724 cn = c * sn
725 if d and _signBit(self.kp2): # PYCHOK no cover
726 cn, dn = dn, cn
727 sn = sn / d # /= chokes PyChecker
728 else:
729 sn = tanh(x)
730 cn = dn = _1_0 / cosh(x)
732 return Elliptic3Tuple(sn, cn, dn, iteration=self._iteration)
734 @Property_RO
735 def _sncndnBulirsch4(self):
736 '''(INTERNAL) Get Bulirsch' 4-tuple C{(c, d, cd, mn_)}.
737 '''
738 # Bulirsch's sncndn routine, p 89.
739 d, mc = 0, self.kp2
740 if _signBit(mc): # PYCHOK no cover
741 d = _1_0 - mc
742 mc = neg(mc / d)
743 d = sqrt(d)
745 mn, a = [], _1_0
746 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC
747 mc = sqrt(mc)
748 mn.append((a, mc))
749 c = (a + mc) * _0_5
750 t = _TolJAC * a
751 if fabs(a - mc) <= t: # 6 trips, quadratic
752 self._iteration += i # accumulate
753 break
754 mc *= a
755 a = c
756 else: # PYCHOK no cover
757 raise _convergenceError(a - mc, t, None, kp=self.kp, kp2=self.kp2)
758 cd = (c * d) if d else c
759 return c, d, cd, tuple(reversed(mn)) # mn reversed!
761 def _sncndnPhi(self, phi):
762 '''(INTERNAL) Helper for C{.fEinv} and C{._fXf}.
763 '''
764 sn, cn = sincos2(phi)
765 return Elliptic3Tuple(sn, cn, self.fDelta(sn, cn))
767 @staticmethod
768 def fRC(x, y):
769 '''Degenerate symmetric integral of the first kind C{RC(x, y)}.
771 @return: C{RC(x, y)}, equivalent to C{RF(x, y, y)}.
773 @see: U{C{RC} definition<https://DLMF.NIST.gov/19.2.E17>} and
774 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
775 '''
776 return _rC(None, x, y)
778 @staticmethod
779 def fRD(x, y, z, *over):
780 '''Degenerate symmetric integral of the third kind C{RD(x, y, z)}.
782 @return: C{RD(x, y, z) / over}, equivalent to C{RJ(x, y, z, z)
783 / over} with C{over} typically 3.
785 @see: U{C{RD} definition<https://DLMF.NIST.gov/19.16.E5>} and
786 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
787 '''
788 return float(_RD(None, x, y, z, *over))
790 @staticmethod
791 def fRF(x, y, z=0):
792 '''Symmetric or complete symmetric integral of the first kind
793 C{RF(x, y, z)} respectively C{RF(x, y)}.
795 @return: C{RF(x, y, z)} or C{RF(x, y)} for missing or zero B{C{z}}.
797 @see: U{C{RF} definition<https://DLMF.NIST.gov/19.16.E1>} and
798 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
799 '''
800 return float(_RF3(None, x, y, z)) if z else _rF2(None, x, y)
802 @staticmethod
803 def fRG(x, y, z=0):
804 '''Symmetric or complete symmetric integral of the second kind
805 C{RG(x, y, z)} respectively C{RG(x, y)}.
807 @return: C{RG(x, y, z)} or C{RG(x, y)} for missing or zero B{C{z}}.
809 @see: U{C{RG} definition<https://DLMF.NIST.gov/19.16.E3>} and
810 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
811 '''
812 G = _RG3(None, x, y, z) if z else (_RG2(None, x, y) * _0_5)
813 return float(G)
815 @staticmethod
816 def fRJ(x, y, z, p):
817 '''Symmetric integral of the third kind C{RJ(x, y, z, p)}.
819 @return: C{RJ(x, y, z, p)}.
821 @see: U{C{RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and
822 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
823 '''
824 return float(_RJ(None, x, y, z, p))
826_allPropertiesOf_n(15, Elliptic) # # PYCHOK assert, see Elliptic.reset
829class EllipticError(_ValueError):
830 '''Elliptic integral, function, convergence or other L{Elliptic} issue.
831 '''
832 pass
835class Elliptic3Tuple(_NamedTuple):
836 '''3-Tuple C{(sn, cn, dn)} all C{scalar}.
837 '''
838 _Names_ = ('sn', 'cn', 'dn')
839 _Units_ = ( Scalar, Scalar, Scalar)
842class _Lxyz(list):
843 '''(INTERNAL) Helper for C{_RD}, C{_RF3} and C{_RJ}.
844 '''
845 _a = None
846 _a0 = None
847# _xyzp = ()
849 def __init__(self, *xyzp): # x, y, z [, p]
850 list.__init__(self, xyzp)
851 self._xyzp = xyzp
853 def a0(self, n):
854 '''Compute the initial C{a}.
855 '''
856 t = tuple(self)
857 m = n - len(t)
858 if m > 0:
859 t += t[-1:] * m
860 try:
861 s = Fsum(*t).fover(n)
862 except ValueError: # Fsum(NAN) exception
863 s = _sum(t) / n
864 self._a0 = self._a = s
865 return s
867 def asr3(self, a):
868 '''Compute next C{a}, C{sqrt(xyz_)} and C{fdot(sqrt(xyz))}.
869 '''
870 L = self
871 # assert a is L._a
872 s = map2(sqrt, L) # sqrt(x), srqt(y), sqrt(z) [, sqrt(p)]
873 try:
874 r = fdot(s[:3], s[1], s[2], s[0]) # sqrt(x) * sqrt(y) + ...
875 except ValueError: # Fsum(NAN) exception
876 r = _sum(s[i] * s[(i + 1) % 3] for i in range(3))
877 L[:] = [(x + r) * _0_25 for x in L]
878 # assert L is self
879 L._a = a = (a + r) * _0_25
880 return a, s, r
882 def Casr4(self, inst, n, Tol, where):
883 '''Yield Carlson 4-tuples C{(m, a, s, r)} plus sentinel,
884 with C{s = sqrt(xyz_)} and C{r = fdot(sqrt(xyz))}.
885 '''
886 L = self
887 a = am = L.a0(n)
888 m = 1
889 q = max(fabs(a - x) for x in L) / Tol
890 for i in range(_TRIPS):
891 if fabs(am) > q: # 5-6 trips
892 _iterations(inst, i)
893 break
894 a, s, r = L.asr3(a)
895 yield m, a, s, r
896 m *= 4
897 am = a * m
898 else: # PYCHOK no cover
899 raise _convergenceError(am, q, where, *L._xyzp, thresh=True)
900 yield m, a, (), None # sentinel: next m, same a, no s, no r
902 def rescale(self, am, *xy_):
903 '''Rescale C{x}, C{y}, ...
904 '''
905 for x in xy_:
906 yield (self._a0 - x) / am
909def _Cab3(inst, x, y, where):
910 '''(INTERNAL) Yield C{(i, a, b)} Carlson 3-tuples.
911 '''
912 a, b = sqrt(x), sqrt(y)
913 if b > a:
914 a, b = b, a
915 yield a, b, 0
916 for i in range(1, _TRIPS):
917 t = _TolRG0 * a
918 if fabs(a - b) <= t: # 3-4 trips
919 _iterations(inst, i - 1)
920 break
921 a, b = ((a + b) * _0_5), sqrt(a * b)
922 yield a, b, i
923 else: # PYCHOK no cover
924 raise _convergenceError(a - b, t, where, x, y)
927def _convergenceError(d, tol, where, *args, **kwds_thresh): # PYCHOK no cover
928 '''(INTERNAL) Return an L{EllipticError}.
929 '''
930 n = Elliptic.__name__
931 if where:
932 n = _DOT_(n, where.__name__)
933 if kwds_thresh:
934 q = _xkwds_pop(kwds_thresh, thresh=False)
935 t = unstr(n, *args, **kwds_thresh)
936 else:
937 q = False
938 t = unstr(n, *args)
939 return EllipticError(Fmt.no_convergence(d, tol, thresh=q), txt=t)
942def _deltaX(sn, cn, dn, cX, _fX):
943 '''(INTERNAL) Helper for C{Elliptic.deltaD} thru C{.deltaPi}.
944 '''
945 if cn is None or dn is None:
946 n = NN(_delta_, _fX.__name__[1:])
947 raise _invokationError(n, sn, cn, dn)
949 if _signBit(cn):
950 cn, sn = -cn, -sn
951 return _fX(sn, cn, dn) * PI_2 / cX - atan2(sn, cn)
954def _Horner(S, e1, E2, E3, E4, E5, *over):
955 '''(INTERNAL) Horner form for C{_RD} and C{_RJ} below.
956 '''
957 E22 = E2**2
958 # Polynomial is <https://DLMF.NIST.gov/19.36.E2>
959 # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52
960 # + 3*E5/26 - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20
961 # + 45*E2**2*E3/272 - 9*(E3*E4+E2*E5)/68)
962 # converted to Horner-like form ...
963 F = Fsum
964 e = e1 * 4084080
965 S *= e
966 S += F(E2 * -540540, 471240).fmul(E5)
967 S += F(E2 * 612612, E3 * -540540, -556920).fmul(E4)
968 S += F(E2 * -706860, E22 * 675675, E3 * 306306, 680680).fmul(E3)
969 S += F(E2 * 417690, E22 * -255255, -875160).fmul(E2)
970 S += 4084080
971 return S.fdiv((over[0] * e) if over else e) # Fsum
974def _invokationError(name, *args): # PYCHOK no cover
975 '''(INTERNAL) Return an L{EllipticError}.
976 '''
977 n = _DOT_(Elliptic.__name__, name)
978 n = _SPACE_(_invokation_, n)
979 return EllipticError(NN(n, repr(args))) # unstr
982def _iterations(inst, i):
983 '''(INTERNAL) Aggregate iterations B{C{i}}.
984 '''
985 if inst and i > 0:
986 inst._iteration += i
989def _3over(a, b):
990 '''(INTERNAL) Return C{3 / (a * b)}.
991 '''
992 return _over(_3_0, a * b)
995def _rC(unused, x, y):
996 '''(INTERNAL) Defined only for y != 0 and x >= 0.
997 '''
998 d = x - y
999 if d < 0: # catch _NaN
1000 # <https://DLMF.NIST.gov/19.2.E18>
1001 d = -d
1002 r = atan(sqrt(d / x)) if x > 0 else PI_2
1003 elif d == _0_0: # XXX d < EPS0? or EPS02 or _EPSmin
1004 d, r = y, _1_0
1005 elif y > 0: # <https://DLMF.NIST.gov/19.2.E19>
1006 r = asinh(sqrt(d / y)) # atanh(sqrt((x - y) / x))
1007 elif y < 0: # <https://DLMF.NIST.gov/19.2.E20>
1008 r = asinh(sqrt(-x / y)) # atanh(sqrt(x / (x - y)))
1009 else:
1010 raise _invokationError(Elliptic.fRC.__name__, x, y)
1011 return r / sqrt(d) # float
1014def _RD(inst, x, y, z, *over):
1015 '''(INTERNAL) Carlson, eqs 2.28 - 2.34.
1016 '''
1017 L = _Lxyz(x, y, z)
1018 S = _Dsum()
1019 for m, a, s, r in L.Casr4(inst, 5, _TolRF, Elliptic.fRD):
1020 if s:
1021 S += _over(_3_0, (z + r) * s[2] * m)
1022 z = L[2] # s[2] = sqrt(z)
1023 x, y = L.rescale(-a * m, x, y)
1024 xy = x * y
1025 z = (x + y) / _3_0
1026 z2 = z**2
1027 return _Horner(S(_1_0), sqrt(a) * a * m,
1028 xy - _6_0 * z2,
1029 (xy * _3_0 - _8_0 * z2) * z,
1030 (xy - z2) * _3_0 * z2,
1031 xy * z2 * z, *over) # Fsum
1034def _rF2(inst, x, y): # 2-arg version, z=0
1035 '''(INTERNAL) Carlson, eqs 2.36 - 2.38.
1036 '''
1037 for a, b, _ in _Cab3(inst, x, y, Elliptic.fRF): # PYCHOK yield
1038 pass
1039 return _over(PI, a + b) # float
1042def _RF3(inst, x, y, z): # 3-arg version
1043 '''(INTERNAL) Carlson, eqs 2.2 - 2.7.
1044 '''
1045 L = _Lxyz(x, y, z)
1046 for m, a, _, _ in L.Casr4(inst, 3, _TolRF, Elliptic.fRF):
1047 pass
1048 x, y = L.rescale(a * m, x, y)
1049 z = neg(x + y)
1050 xy = x * y
1051 e2 = xy - z**2
1052 e3 = xy * z
1053 e4 = e2**2
1054 # Polynomial is <https://DLMF.NIST.gov/19.36.E1>
1055 # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44
1056 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16)
1057 # converted to Horner-like form ...
1058 S = Fsum(e4 * 15015, e3 * 6930, e2 * -16380, 17160).fmul(e3)
1059 S += Fsum(e4 * -5775, e2 * 10010, -24024).fmul(e2)
1060 S += 240240
1061 return S.fdiv(sqrt(a) * 240240) # Fsum
1064def _RG2(inst, x, y): # 2-args and I{doubled}
1065 '''(INTERNAL) Carlson, eqs 2.36 - 2.39.
1066 '''
1067 m = -1 # neg!
1068 S = _Dsum()
1069 for a, b, i in _Cab3(inst, x, y, Elliptic.fRG):
1070 if i:
1071 S += (a - b)**2 * m
1072 m *= 2
1073 else:
1074 S += (a + b)**2 * _0_5
1075 return S(PI_2 / (a + b)) # Fsum
1078def _RG3(inst, x, y, z): # 3-arg version
1079 '''(INTERNAL) Never called with zero B{C{z}}, see C{.fRG}.
1080 '''
1081# if not z:
1082# y, z = z, y
1083 R = _RF3(inst, x, y, z)
1084 rd = (x - z) * (z - y) # - (y - z)
1085 if rd: # Carlson, eq 1.7
1086 R += _RD(inst, x, y, z, _3_0 * z / rd)
1087 r = x * y
1088 if r:
1089 R += sqrt(r / z**3)
1090 return R.fmul(z * _0_5) # Fsum
1093def _RJ(inst, x, y, z, p, *over):
1094 '''(INTERNAL) Carlson, eqs 2.17 - 2.25.
1095 '''
1096 def _xyzp(x, y, z, p):
1097 return (x + p) * (y + p) * (z + p)
1099 L = _Lxyz(x, y, z, p)
1100 n = neg(_xyzp(x, y, z, -p))
1101 S = _Dsum()
1102 for m, a, s, _ in L.Casr4(inst, 5, _TolRD, Elliptic.fRJ):
1103 if s:
1104 d = _xyzp(*s)
1105 if n:
1106 rc = _rC(inst, _1_0, n / d**2 + _1_0)
1107 n *= _1_64th # /= chokes PyChecker
1108 else:
1109 rc = _1_0 # == _rC(None, _1_0, _1_0)
1110 S += rc / (d * m)
1111 x, y, z = L.rescale(a * m, x, y, z)
1112 xyz = x * y * z
1113 p = Fsum(x, y, z).fover(_N_2_0)
1114 p2 = p**2
1115 p3 = p2 * p
1116 E2 = Fsum(x * y, x * z, y * z, -p2 * _3_0)
1117 E2p = E2 * p
1118 return _Horner(S(_6_0), sqrt(a) * a * m, E2,
1119 Fsum(p3 * _4_0, xyz, E2p * _2_0),
1120 Fsum(p3 * _3_0, E2p, xyz * _2_0).fmul(p),
1121 xyz * p2, *over) # Fsum
1123# **) MIT License
1124#
1125# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1126#
1127# Permission is hereby granted, free of charge, to any person obtaining a
1128# copy of this software and associated documentation files (the "Software"),
1129# to deal in the Software without restriction, including without limitation
1130# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1131# and/or sell copies of the Software, and to permit persons to whom the
1132# Software is furnished to do so, subject to the following conditions:
1133#
1134# The above copyright notice and this permission notice shall be included
1135# in all copies or substantial portions of the Software.
1136#
1137# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1138# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1139# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1140# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1141# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1142# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1143# OTHER DEALINGS IN THE SOFTWARE.