Coverage for pygeodesy/elliptic.py: 99%
455 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-06 15:27 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-06 15:27 -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:Charles@Karney.com>} (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, _EPStol as _TolJAC, INF, PI, PI_2, PI_4, \
79 _0_0, _0_125, _0_25, _0_5, _1_0, _1_64th, _2_0, \
80 _N_2_0, _3_0, _4_0, _6_0, _8_0, _180_0, _360_0
81from pygeodesy.errors import _ValueError, _xkwds_pop
82from pygeodesy.fmath import fdot, hypot1, Fsum
83# from pygeodesy.fsums import Fsum # from .fmath
84from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _SPACE_
85from pygeodesy.karney import _ALL_LAZY, _signBit
86# from pygeodesy.lazily import _ALL_LAZY # from .karney
87from pygeodesy.named import _Named, _NamedTuple, _NotImplemented
88from pygeodesy.props import _allPropertiesOf_n, Property_RO, property_RO, \
89 _update_all
90from pygeodesy.streprs import Fmt, unstr
91from pygeodesy.units import Scalar, Scalar_
92from pygeodesy.utily import sincos2, sincos2d
94from math import asinh, atan, atan2, ceil, cosh, fabs, floor, sin, sqrt, tanh
96__all__ = _ALL_LAZY.elliptic
97__version__ = '23.07.07'
99_invokation_ = 'invokation'
100_TolRD = pow(EPS * 0.002, _0_125) # 8th root: quadquadratic?, octic?, ocrt?
101_TolRF = pow(EPS * 0.030, _0_125) # 4th root: biquadratic, quartic, qurt?
102_TolRG0 = _TolJAC * 2.7
103_TRIPS = 31 # Max depth, 7 might be sufficient
106class _Complete(object):
107 '''(INTERAL) Hold complete integrals.
108 '''
109 def __init__(self, **kwds):
110 self.__dict__ = kwds
113class _Deferred(list):
114 '''(INTERNAL) Collector for L{_Deferred_Fsum}.
115 '''
116 def __init__(self, *xs):
117 if xs:
118 list.__init__(self, xs)
120 def __add__(self, other): # PYCHOK no cover
121 return _NotImplemented(self, other)
123 def __iadd__(self, x): # overide C{list} += C{list}
124 # assert isscalar(x)
125 self.append(float(x))
126 return self
128 def __imul__(self, other): # PYCHOK no cover
129 return _NotImplemented(self, other)
131 def __isub__(self, x): # PYCHOK no cover
132 # assert isscalar(x)
133 self.append(-float(x))
134 return self
136 def __mul__(self, other): # PYCHOK no cover
137 return _NotImplemented(self, other)
139 def __radd__(self, other): # PYCHOK no cover
140 return _NotImplemented(self, other)
142 def __rsub__(self, other): # PYCHOK no cover
143 return _NotImplemented(self, other)
145 def __sub__(self, other): # PYCHOK no cover
146 return _NotImplemented(self, other)
148 @property_RO
149 def Fsum(self):
150 # get a L{_Deferred_Fsum} instance, pre-named
151 return _Deferred_Fsum()._facc(self) # known C{float}s
154class _Deferred_Fsum(Fsum):
155 '''(INTERNAL) Deferred L{Fsum}.
156 '''
157 name = NN # pre-named, overridden below
159 def _update(self, **other): # PYCHOK don't ...
160 # ... waste time zapping non-existing Property/_ROs
161 if other or len(self.__dict__) > 2:
162 Fsum._update(self, **other)
164_Deferred_Fsum.name = _Deferred_Fsum.__name__ # PYCHOK once
167class Elliptic(_Named):
168 '''Elliptic integrals and functions.
170 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/
171 html/classGeographicLib_1_1EllipticFunction.html#details>}.
172 '''
173 _alpha2 = 0
174 _alphap2 = 0
175 _eps = EPS
176 _k2 = 0
177 _kp2 = 0
179 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None, name=NN):
180 '''Constructor, specifying the C{modulus} and C{parameter}.
182 @kwarg name: Optional name (C{str}).
184 @see: Method L{Elliptic.reset} for further details.
186 @note: If only elliptic integrals of the first and second kinds
187 are needed, use C{B{alpha2}=0}, the default value. In
188 that case, we have C{Π(φ, 0, k) = F(φ, k), G(φ, 0, k) =
189 E(φ, k)} and C{H(φ, 0, k) = F(φ, k) - D(φ, k)}.
190 '''
191 self.reset(k2=k2, alpha2=alpha2, kp2=kp2, alphap2=alphap2)
193 if name:
194 self.name = name
196 @Property_RO
197 def alpha2(self):
198 '''Get α^2, the square of the parameter (C{float}).
199 '''
200 return self._alpha2
202 @Property_RO
203 def alphap2(self):
204 '''Get α'^2, the square of the complementary parameter (C{float}).
205 '''
206 return self._alphap2
208 @Property_RO
209 def cD(self):
210 '''Get Jahnke's complete integral C{D(k)} (C{float}),
211 U{defined<https://DLMF.NIST.gov/19.2.E6>}.
212 '''
213 return self._reset_cDcEcKcKE_eps.cD
215 @Property_RO
216 def cE(self):
217 '''Get the complete integral of the second kind C{E(k)}
218 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}.
219 '''
220 return self._reset_cDcEcKcKE_eps.cE
222 @Property_RO
223 def cG(self):
224 '''Get Legendre's complete geodesic longitude integral
225 C{G(α^2, k)} (C{float}).
226 '''
227 return self._reset_cGcHcPi.cG
229 @Property_RO
230 def cH(self):
231 '''Get Cayley's complete geodesic longitude difference integral
232 C{H(α^2, k)} (C{float}).
233 '''
234 return self._reset_cGcHcPi.cH
236 @Property_RO
237 def cK(self):
238 '''Get the complete integral of the first kind C{K(k)}
239 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}.
240 '''
241 return self._reset_cDcEcKcKE_eps.cK
243 @Property_RO
244 def cKE(self):
245 '''Get the difference between the complete integrals of the
246 first and second kinds, C{K(k) − E(k)} (C{float}).
247 '''
248 return self._reset_cDcEcKcKE_eps.cKE
250 @Property_RO
251 def cPi(self):
252 '''Get the complete integral of the third kind C{Pi(α^2, k)}
253 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}.
254 '''
255 return self._reset_cGcHcPi.cPi
257 def deltaD(self, sn, cn, dn):
258 '''The periodic Jahnke's incomplete elliptic integral.
260 @arg sn: sin(φ).
261 @arg cn: cos(φ).
262 @arg dn: sqrt(1 − k2 * sin(2φ)).
264 @return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}).
266 @raise EllipticError: Invalid invokation or no convergence.
267 '''
268 return self._deltaX(sn, cn, dn, self.cD, self.fD)
270 def deltaE(self, sn, cn, dn):
271 '''The periodic incomplete integral of the second kind.
273 @arg sn: sin(φ).
274 @arg cn: cos(φ).
275 @arg dn: sqrt(1 − k2 * sin(2φ)).
277 @return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}).
279 @raise EllipticError: Invalid invokation or no convergence.
280 '''
281 return self._deltaX(sn, cn, dn, self.cE, self.fE)
283 def deltaEinv(self, stau, ctau):
284 '''The periodic inverse of the incomplete integral of the second kind.
286 @arg stau: sin(τ)
287 @arg ctau: cos(τ)
289 @return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}).
291 @raise EllipticError: No convergence.
292 '''
293 # Function is periodic with period pi
294 t = atan2(-stau, -ctau) if _signBit(ctau) else atan2(stau, ctau)
295 return self.fEinv(t * self.cE / PI_2) - t
297 def deltaF(self, sn, cn, dn):
298 '''The periodic incomplete integral of the first kind.
300 @arg sn: sin(φ).
301 @arg cn: cos(φ).
302 @arg dn: sqrt(1 − k2 * sin(2φ)).
304 @return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}).
306 @raise EllipticError: Invalid invokation or no convergence.
307 '''
308 return self._deltaX(sn, cn, dn, self.cK, self.fF)
310 def deltaG(self, sn, cn, dn):
311 '''Legendre's periodic geodesic longitude integral.
313 @arg sn: sin(φ).
314 @arg cn: cos(φ).
315 @arg dn: sqrt(1 − k2 * sin(2φ)).
317 @return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}).
319 @raise EllipticError: Invalid invokation or no convergence.
320 '''
321 return self._deltaX(sn, cn, dn, self.cG, self.fG)
323 def deltaH(self, sn, cn, dn):
324 '''Cayley's periodic geodesic longitude difference integral.
326 @arg sn: sin(φ).
327 @arg cn: cos(φ).
328 @arg dn: sqrt(1 − k2 * sin(2φ)).
330 @return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}).
332 @raise EllipticError: Invalid invokation or no convergence.
333 '''
334 return self._deltaX(sn, cn, dn, self.cH, self.fH)
336 def deltaPi(self, sn, cn, dn):
337 '''The periodic incomplete integral of the third kind.
339 @arg sn: sin(φ).
340 @arg cn: cos(φ).
341 @arg dn: sqrt(1 − k2 * sin(2φ)).
343 @return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ
344 (C{float}).
346 @raise EllipticError: Invalid invokation or no convergence.
347 '''
348 return self._deltaX(sn, cn, dn, self.cPi, self.fPi)
350 def _deltaX(self, sn, cn, dn, cX, fX):
351 '''(INTERNAL) Helper for C{.deltaD} thru C{.deltaPi}.
352 '''
353 if cn is None or dn is None:
354 n = NN(_delta_, fX.__name__[1:])
355 raise _invokationError(n, sn, cn, dn)
357 if _signBit(cn):
358 cn, sn = -cn, -sn
359 return fX(sn, cn, dn) * PI_2 / cX - atan2(sn, cn)
361 @Property_RO
362 def eps(self):
363 '''Get epsilon (C{float}).
364 '''
365 return self._reset_cDcEcKcKE_eps.eps
367 def fD(self, phi_or_sn, cn=None, dn=None):
368 '''Jahnke's incomplete elliptic integral in terms of
369 Jacobi elliptic functions.
371 @arg phi_or_sn: φ or sin(φ).
372 @kwarg cn: C{None} or cos(φ).
373 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
375 @return: D(φ, k) as though φ ∈ (−π, π] (C{float}),
376 U{defined<https://DLMF.NIST.gov/19.2.E4>}.
378 @raise EllipticError: Invalid invokation or no convergence.
379 '''
380 def _fD(sn, cn, dn):
381 r = fabs(sn)**3
382 if r:
383 r = _RD(self, cn**2, dn**2, _1_0, _3_0 / r)
384 return r
386 return self._fXf(phi_or_sn, cn, dn, self.cD,
387 self.deltaD, _fD)
389 def fDelta(self, sn, cn):
390 '''The C{Delta} amplitude function.
392 @arg sn: sin(φ).
393 @arg cn: cos(φ).
395 @return: sqrt(1 − k2 * sin(2φ)) (C{float}).
396 '''
397 k2 = self.k2
398 s = (_1_0 - k2 * sn**2) if k2 < 0 else (self.kp2
399 + ((k2 * cn**2) if k2 > 0 else _0_0))
400 return sqrt(s) if s else _0_0
402 def fE(self, phi_or_sn, cn=None, dn=None):
403 '''The incomplete integral of the second kind in terms of
404 Jacobi elliptic functions.
406 @arg phi_or_sn: φ or sin(φ).
407 @kwarg cn: C{None} or cos(φ).
408 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
410 @return: E(φ, k) as though φ ∈ (−π, π] (C{float}),
411 U{defined<https://DLMF.NIST.gov/19.2.E5>}.
413 @raise EllipticError: Invalid invokation or no convergence.
414 '''
415 def _fE(sn, cn, dn):
416 '''(INTERNAL) Core of C{.fE}.
417 '''
418 if sn:
419 sn2, cn2, dn2 = sn**2, cn**2, dn**2
420 kp2, k2 = self.kp2, self.k2
421 if k2 <= 0: # Carlson, eq. 4.6, <https://DLMF.NIST.gov/19.25.E9>
422 ei = _RF3(self, cn2, dn2, _1_0)
423 if k2:
424 ei -= _RD(self, cn2, dn2, _1_0, _3over(k2, sn2))
425 elif kp2 >= 0: # <https://DLMF.NIST.gov/19.25.E10>
426 ei = k2 * fabs(cn) / dn
427 if kp2:
428 ei += (_RD( self, cn2, _1_0, dn2, _3over(k2, sn2)) +
429 _RF3(self, cn2, dn2, _1_0)) * kp2
430 else: # <https://DLMF.NIST.gov/19.25.E11>
431 ei = dn / fabs(cn)
432 ei -= _RD(self, dn2, _1_0, cn2, _3over(kp2, sn2))
433 ei *= fabs(sn)
434 else: # PYCHOK no cover
435 ei = _0_0
436 return ei
438 return self._fXf(phi_or_sn, cn, dn, self.cE,
439 self.deltaE, _fE)
441 def fEd(self, deg):
442 '''The incomplete integral of the second kind with
443 the argument given in degrees.
445 @arg deg: Angle (C{degrees}).
447 @return: E(π B{C{deg}}/180, k) (C{float}).
449 @raise EllipticError: No convergence.
450 '''
451 if fabs(deg) < _180_0:
452 e = _0_0
453 else: # PYCHOK no cover
454 e = ceil(deg / _360_0 - _0_5)
455 deg -= e * _360_0
456 e *= self.cE * _4_0
457 sn, cn = sincos2d(deg)
458 return self.fE(sn, cn, self.fDelta(sn, cn)) + e
460 def fEinv(self, x):
461 '''The inverse of the incomplete integral of the second kind.
463 @arg x: Argument (C{float}).
465 @return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}}
466 (C{float}).
468 @raise EllipticError: No convergence.
469 '''
470 E2 = self.cE * _2_0
471 n = floor(x / E2 + _0_5)
472 y = x - E2 * n # y now in [-ec, ec)
473 # linear approximation
474 phi = PI * y / E2 # phi in [-pi/2, pi/2)
475 Phi = Fsum(phi)
476 # first order correction
477 phi = Phi.fsum_(self.eps * sin(phi * _2_0) / _N_2_0)
478 # For kp2 close to zero use asin(x/.cE) or J. P. Boyd,
479 # Applied Math. and Computation 218, 7005-7013 (2012)
480 # <https://DOI.org/10.1016/j.amc.2011.12.021>
481 fE = self.fE
482 _sncndnPhi = self._sncndnPhi
483 _Phi2 = Phi.fsum2_
484 self._iteration = 0 # aggregate
485 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC
486 sn, cn, dn = _sncndnPhi(phi)
487 phi, e = _Phi2((y - fE(sn, cn, dn)) / dn)
488 if fabs(e) < _TolJAC:
489 self._iteration += i
490 break
491 else: # PYCHOK no cover
492 raise _no_convergenceError(e, _TolJAC, self.fEinv, x)
493 return Phi.fsum_(n * PI) if n else phi
495 def fF(self, phi_or_sn, cn=None, dn=None):
496 '''The incomplete integral of the first kind in terms of
497 Jacobi elliptic functions.
499 @arg phi_or_sn: φ or sin(φ).
500 @kwarg cn: C{None} or cos(φ).
501 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
503 @return: F(φ, k) as though φ ∈ (−π, π] (C{float}),
504 U{defined<https://DLMF.NIST.gov/19.2.E4>}.
506 @raise EllipticError: Invalid invokation or no convergence.
507 '''
508 def _fF(sn, cn, dn):
509 r = fabs(sn)
510 if r:
511 r *= _RF3(self, cn**2, dn**2, _1_0)
512 return r
514 return self._fXf(phi_or_sn, cn, dn, self.cK,
515 self.deltaF, _fF)
517 def fG(self, phi_or_sn, cn=None, dn=None):
518 '''Legendre's geodesic longitude integral in terms of
519 Jacobi elliptic functions.
521 @arg phi_or_sn: φ or sin(φ).
522 @kwarg cn: C{None} or cos(φ).
523 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
525 @return: G(φ, k) as though φ ∈ (−π, π] (C{float}).
527 @raise EllipticError: Invalid invokation or no convergence.
529 @note: Legendre expresses the longitude of a point on the
530 geodesic in terms of this combination of elliptic
531 integrals in U{Exercices de Calcul Intégral, Vol 1
532 (1811), p 181<https://Books.Google.com/books?id=
533 riIOAAAAQAAJ&pg=PA181>}.
535 @see: U{Geodesics in terms of elliptic integrals<https://
536 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>}
537 for the expression for the longitude in terms of this function.
538 '''
539 return self._fXa(phi_or_sn, cn, dn, self.alpha2 - self.k2,
540 self.cG, self.deltaG)
542 def fH(self, phi_or_sn, cn=None, dn=None):
543 '''Cayley's geodesic longitude difference integral in terms of
544 Jacobi elliptic functions.
546 @arg phi_or_sn: φ or sin(φ).
547 @kwarg cn: C{None} or cos(φ).
548 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
550 @return: H(φ, k) as though φ ∈ (−π, π] (C{float}).
552 @raise EllipticError: Invalid invokation or no convergence.
554 @note: Cayley expresses the longitude difference of a point
555 on the geodesic in terms of this combination of
556 elliptic integrals in U{Phil. Mag. B{40} (1870), p 333
557 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}.
559 @see: U{Geodesics in terms of elliptic integrals<https://
560 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>}
561 for the expression for the longitude in terms of this function.
562 '''
563 return self._fXa(phi_or_sn, cn, dn, -self.alphap2,
564 self.cH, self.deltaH)
566 def fPi(self, phi_or_sn, cn=None, dn=None):
567 '''The incomplete integral of the third kind in terms of
568 Jacobi elliptic functions.
570 @arg phi_or_sn: φ or sin(φ).
571 @kwarg cn: C{None} or cos(φ).
572 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
574 @return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}).
576 @raise EllipticError: Invalid invokation or no convergence.
577 '''
578 if dn is None and cn is not None: # and isscalar(phi_or_sn)
579 dn = self.fDelta(phi_or_sn, cn) # in .triaxial
580 return self._fXa(phi_or_sn, cn, dn, self.alpha2,
581 self.cPi, self.deltaPi)
583 def _fXa(self, phi_or_sn, cn, dn, aX, cX, deltaX):
584 '''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}.
585 '''
586 def _fX(sn, cn, dn):
587 if sn:
588 cn2, sn2, dn2 = cn**2, sn**2, dn**2
589 r = _RF3(self, cn2, dn2, _1_0)
590 if aX:
591 z = cn2 + sn2 * self.alphap2
592 r += _RJ(self, cn2, dn2, _1_0, z, _3over(aX, sn2))
593 r *= fabs(sn)
594 else: # PYCHOK no cover
595 r = _0_0
596 return r
598 return self._fXf(phi_or_sn, cn, dn, cX, deltaX, _fX)
600 def _fXf(self, phi_or_sn, cn, dn, cX, deltaX, fX):
601 '''(INTERNAL) Helper for C{f.D}, C{.fE}, C{.fF} and C{._fXa}.
602 '''
603 self._iteration = 0 # aggregate
604 phi = sn = phi_or_sn
605 if cn is dn is None: # fX(phi) call
606 sn, cn, dn = self._sncndnPhi(phi)
607 if fabs(phi) >= PI: # PYCHOK no cover
608 return (deltaX(sn, cn, dn) + phi) * cX / PI_2
609 # fall through
610 elif cn is None or dn is None:
611 n = NN(_f_, deltaX.__name__[5:])
612 raise _invokationError(n, sn, cn, dn)
614 if _signBit(cn): # enforce usual trig-like symmetries
615 xi = _2_0 * cX - fX(sn, cn, dn)
616 elif cn > 0:
617 xi = fX(sn, cn, dn)
618 else:
619 xi = cX
620 return copysign0(xi, sn)
622 @Property_RO
623 def k2(self):
624 '''Get k^2, the square of the modulus (C{float}).
625 '''
626 return self._k2
628 @Property_RO
629 def kp2(self):
630 '''Get k'^2, the square of the complementary modulus (C{float}).
631 '''
632 return self._kp2
634 def reset(self, k2=0, alpha2=0, kp2=None, alphap2=None): # MCCABE 13
635 '''Reset the modulus, parameter and the complementaries.
637 @kwarg k2: Modulus squared (C{float}, NINF <= k^2 <= 1).
638 @kwarg alpha2: Parameter squared (C{float}, NINF <= α^2 <= 1).
639 @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0).
640 @kwarg alphap2: Complementary parameter squared (C{float}, α'^2 >= 0).
642 @raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}}
643 or B{C{alphap2}}.
645 @note: The arguments must satisfy C{B{k2} + B{kp2} = 1} and
646 C{B{alpha2} + B{alphap2} = 1}. No checking is done
647 that these conditions are met to enable accuracy to be
648 maintained, e.g., when C{k} is very close to unity.
649 '''
650 _update_all(self, _Named.iteration._uname, Base=Property_RO)
652 self._k2 = Scalar_(k2=k2, Error=EllipticError, low=None, high=_1_0)
653 self._kp2 = Scalar_(kp2=((_1_0 - k2) if kp2 is None else kp2), Error=EllipticError)
655 self._alpha2 = Scalar_(alpha2=alpha2, Error=EllipticError, low=None, high=_1_0)
656 self._alphap2 = Scalar_(alphap2=((_1_0 - alpha2) if alphap2 is None else alphap2),
657 Error=EllipticError)
659 # Values of complete elliptic integrals for k = 0,1 and alpha = 0,1
660 # K E D
661 # k = 0: pi/2 pi/2 pi/4
662 # k = 1: inf 1 inf
663 # Pi G H
664 # k = 0, alpha = 0: pi/2 pi/2 pi/4
665 # k = 1, alpha = 0: inf 1 1
666 # k = 0, alpha = 1: inf inf pi/2
667 # k = 1, alpha = 1: inf inf inf
668 #
669 # G(0, k) = Pi(0, k) = H(1, k) = E(k)
670 # H(0, k) = K(k) - D(k)
671 # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2))
672 # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1))
673 # Pi(alpha2, 1) = inf
674 # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2)
676 @Property_RO
677 def _reset_cDcEcKcKE_eps(self):
678 '''(INTERNAL) Get the complete integrals D, E, K and KE plus C{eps}.
679 '''
680 k2 = self.k2
681 if k2:
682 kp2 = self.kp2
683 if kp2:
684 self._iteration = 0
685 # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3
686 # <https://DLMF.NIST.gov/19.25.E1>
687 cD = _RD(self, _0_0, kp2, _1_0, _3_0)
688 # Complete elliptic integral E(k), Carlson eq. 4.2
689 # <https://DLMF.NIST.gov/19.25.E1>
690 cE = _RG2(self, kp2, _1_0)
691 # Complete elliptic integral K(k), Carlson eq. 4.1
692 # <https://DLMF.NIST.gov/19.25.E1>
693 cK = _RF2(self, kp2, _1_0)
694 cKE = k2 * cD
695 eps = k2 / (sqrt(kp2) + _1_0)**2
696 else: # PYCHOK no cover
697 cD = cK = cKE = INF
698 cE = _1_0
699 eps = k2
700 else: # PYCHOK no cover
701 cD = PI_4
702 cE = cK = PI_2
703 cKE = _0_0 # k2 * cD
704 eps = EPS
706 return _Complete(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps)
708 @Property_RO
709 def _reset_cGcHcPi(self):
710 '''(INTERNAL) Get the complete integrals G, H and Pi.
711 '''
712 self._iteration = 0
713 alpha2 = self.alpha2
714 if alpha2:
715 alphap2 = self.alphap2
716 if alphap2:
717 kp2 = self.kp2
718 if kp2: # <https://DLMF.NIST.gov/19.25.E2>
719 rj = _RJ(self, _0_0, kp2, _1_0, alphap2, _3_0)
720 cPi = cH = cG = self.cK
721 cG += (alpha2 - self.k2) * rj # G(alpha2, k)
722 cH -= alphap2 * rj # H(alpha2, k)
723 cPi += alpha2 * rj # Pi(alpha2, k)
724 else: # PYCHOK no cover
725 cG = cH = _RC(self, _1_0, alphap2)
726 cPi = INF # XXX or NAN?
727 else: # PYCHOK no cover
728 cG = cH = cPi = INF # XXX or NAN?
729 else:
730 cG, cPi, kp2 = self.cE, self.cK, self.kp2
731 # H = K - D but this involves large cancellations if k2 is near 1.
732 # So write (for alpha2 = 0)
733 # H = int(cos(phi)**2/sqrt(1-k2*sin(phi)**2),phi,0,pi/2)
734 # = 1/sqrt(1-k2) * int(sin(phi)**2/sqrt(1-k2/kp2*sin(phi)**2,...)
735 # = 1/kp * D(i*k/kp)
736 # and use D(k) = RD(0, kp2, 1) / 3
737 # so H = 1/kp * RD(0, 1/kp2, 1) / 3
738 # = kp2 * RD(0, 1, kp2) / 3
739 # using <https://DLMF.NIST.gov/19.20.E18>. Equivalently
740 # RF(x, 1) - RD(0, x, 1)/3 = x * RD(0, 1, x)/3 for x > 0
741 # For k2 = 1 and alpha2 = 0, we have
742 # H = int(cos(phi),...) = 1
743 cH = _RD(self, _0_0, _1_0, kp2, _3_0 / kp2) if kp2 else _1_0
745 return _Complete(cG=cG, cH=cH, cPi=cPi)
747 def sncndn(self, x):
748 '''The Jacobi elliptic function.
750 @arg x: The argument (C{float}).
752 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with
753 C{*n(B{x}, k)}.
755 @raise EllipticError: No convergence.
756 '''
757 self._iteration = 0 # reset
758 # Bulirsch's sncndn routine, p 89.
759 if self.kp2:
760 c, d, cd, mn_ = self._sncndnBulirsch4
761 dn = _1_0
762 sn, cn = sincos2(x * cd)
763 if sn:
764 a = cn / sn
765 c *= a
766 for m, n in mn_:
767 a *= c
768 c *= dn
769 dn = (n + a) / (m + a)
770 a = c / m
771 sn = copysign0(_1_0 / hypot1(c), sn) # _signBit(sn)
772 cn = c * sn
773 if d and _signBit(self.kp2): # PYCHOK no cover
774 cn, dn = dn, cn
775 sn = sn / d # /= chokes PyChecker
776 else:
777 sn = tanh(x)
778 cn = dn = _1_0 / cosh(x)
780 return Elliptic3Tuple(sn, cn, dn, iteration=self._iteration)
782 @Property_RO
783 def _sncndnBulirsch4(self):
784 '''(INTERNAL) Get Bulirsch' 4-tuple C{(c, d, cd, mn_)}.
785 '''
786 # Bulirsch's sncndn routine, p 89.
787 d, mc = 0, self.kp2
788 if _signBit(mc): # PYCHOK no cover
789 d = _1_0 - mc
790 mc = neg(mc / d)
791 d = sqrt(d)
793 mn, a = [], _1_0
794 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC
795 # This converges quadratically, max 6 trips
796 mc = sqrt(mc)
797 mn.append((a, mc))
798 c = (a + mc) * _0_5
799 t = _TolJAC * a
800 if fabs(a - mc) <= t:
801 self._iteration += i # accumulate
802 break
803 mc *= a
804 a = c
805 else: # PYCHOK no cover
806 raise _no_convergenceError(a - mc, t, None, kp=self.kp, kp2=self.kp2)
807 cd = (c * d) if d else c
808 return c, d, cd, tuple(reversed(mn)) # mn reversed!
810 def _sncndnPhi(self, phi):
811 '''(INTERNAL) Helper for C{.fEinv} and C{._fXf}.
812 '''
813 sn, cn = sincos2(phi)
814 return Elliptic3Tuple(sn, cn, self.fDelta(sn, cn))
816 @staticmethod
817 def fRC(x, y):
818 '''Degenerate symmetric integral of the first kind C{RC(x, y)}.
820 @return: C{RC(x, y)}, equivalent to C{RF(x, y, y)}.
822 @see: U{C{RC} definition<https://DLMF.NIST.gov/19.2.E17>} and
823 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
824 '''
825 return _RC(None, x, y)
827 @staticmethod
828 def fRD(x, y, z, *over):
829 '''Degenerate symmetric integral of the third kind C{RD(x, y, z)}.
831 @return: C{RD(x, y, z) / over}, equivalent to C{RJ(x, y, z, z)
832 / over} with C{over} typically 3.
834 @see: U{C{RD} definition<https://DLMF.NIST.gov/19.16.E5>} and
835 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
836 '''
837 return _RD(None, x, y, z, *over)
839 @staticmethod
840 def fRF(x, y, z=0):
841 '''Symmetric or complete symmetric integral of the first kind
842 C{RF(x, y, z)} respectively C{RF(x, y)}.
844 @return: C{RF(x, y, z)} or C{RF(x, y)} for missing or zero B{C{z}}.
846 @see: U{C{RF} definition<https://DLMF.NIST.gov/19.16.E1>} and
847 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
848 '''
849 return _RF3(None, x, y, z) if z else _RF2(None, x, y)
851 @staticmethod
852 def fRG(x, y, z=0):
853 '''Symmetric or complete symmetric integral of the second kind
854 C{RG(x, y, z)} respectively C{RG(x, y)}.
856 @return: C{RG(x, y, z)} or C{RG(x, y)} for missing or zero B{C{z}}.
858 @see: U{C{RG} definition<https://DLMF.NIST.gov/19.16.E3>} and
859 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
860 '''
861 return _RG3(None, x, y, z) if z else (_RG2(None, x, y) * _0_5)
863 @staticmethod
864 def fRJ(x, y, z, p):
865 '''Symmetric integral of the third kind C{RJ(x, y, z, p)}.
867 @return: C{RJ(x, y, z, p)}.
869 @see: U{C{RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and
870 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
871 '''
872 return _RJ(None, x, y, z, p)
874_allPropertiesOf_n(15, Elliptic) # # PYCHOK assert, see Elliptic.reset
877class EllipticError(_ValueError):
878 '''Elliptic integral, function, convergence or other L{Elliptic} issue.
879 '''
880 pass
883class Elliptic3Tuple(_NamedTuple):
884 '''3-Tuple C{(sn, cn, dn)} all C{scalar}.
885 '''
886 _Names_ = ('sn', 'cn', 'dn')
887 _Units_ = ( Scalar, Scalar, Scalar)
890class _Lxyz(list):
891 '''(INTERNAL) Helper for C{_RD}, C{_RF3} and C{_RJ}.
892 '''
893 _a = None
894 _a0 = None
896 def __init__(self, *xyz_): # x, y, z [, p]
897 list.__init__(self, xyz_)
899 def a0(self, n):
900 '''Compute the initial C{a}.
901 '''
902 t = tuple(self)
903 m = n - len(t)
904 if m > 0:
905 t += t[-1:] * m
906 try:
907 s = Fsum(*t).fover(n)
908 except ValueError: # NAN
909 s = sum(t) / n
910 self._a0 = self._a = s
911 return s
913 def asr3(self, a):
914 '''Compute next C{a}, C{sqrt(xyz_)} and C{fdot(sqrt(xyz))}.
915 '''
916 L = self
917 # assert a is L._a
918 s = map2(sqrt, L) # sqrt(x), srqt(y), sqrt(z) [, sqrt(p)]
919 try:
920 r = fdot(s[:3], s[1], s[2], s[0]) # sqrt(x) * sqrt(y) + ...
921 except ValueError: # NAN
922 r = sum(s[i] * s[(i + 1) % 3] for i in range(3))
923 L[:] = [(x + r) * _0_25 for x in L]
924 # assert L is self
925 L._a = a = (a + r) * _0_25
926 return a, s, r
928 def rescale(self, am, *xy_):
929 '''Rescale C{x}, C{y}, ...
930 '''
931 for x in xy_:
932 yield (self._a0 - x) / am
934 def thresh(self, Tol):
935 '''Return the convergence threshold.
936 '''
937 return max(fabs(self._a0 - x) for x in self) / Tol
940def _horner(S, e1, E2, E3, E4, E5, *over):
941 '''(INTERNAL) Horner form for C{_RD} and C{_RJ} below.
942 '''
943 E22 = E2**2
944 # Polynomial is <https://DLMF.NIST.gov/19.36.E2>
945 # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52
946 # + 3*E5/26 - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20
947 # + 45*E2**2*E3/272 - 9*(E3*E4+E2*E5)/68)
948 # converted to Horner-like form ...
949 e = e1 * 4084080
950 S *= e
951 S += Fsum( E2 * -540540, 471240).fmul(E5)
952 S += Fsum( E3 * -540540, E2 * 612612, -556920).fmul(E4)
953 S += Fsum(E22 * 675675, E3 * 306306, E2 * -706860, 680680).fmul(E3)
954 S += Fsum(E22 * -255255, E2 * 417690, -875160).fmul(E2)
955 S += 4084080
956 return S.fover((e * over[0]) if over else e)
959def _invokationError(name, *args): # PYCHOK no cover
960 '''(INTERNAL) Return an L{EllipticError}.
961 '''
962 n = _DOT_(Elliptic.__name__, name)
963 n = _SPACE_(_invokation_, n)
964 return EllipticError(NN(n, repr(args))) # unstr
967def _iterations(inst, i):
968 '''(INTERNAL) Aggregate iterations B{C{i}}.
969 '''
970 if inst:
971 inst._iteration += i
974def _no_convergenceError(d, tol, where, *args, **kwds_thresh): # PYCHOK no cover
975 '''(INTERNAL) Return an L{EllipticError}.
976 '''
977 n = Elliptic.__name__
978 if where:
979 n = _DOT_(n, where.__name__)
980 if kwds_thresh:
981 q = _xkwds_pop(kwds_thresh, thresh=False)
982 t = unstr(n, *args, **kwds_thresh)
983 else:
984 q = False
985 t = unstr(n, *args)
986 return EllipticError(Fmt.no_convergence(d, tol, thresh=q), txt=t)
989def _3over(a, b):
990 '''(INTERNAL) Return C{3 / (a * b)}.
991 '''
992 return _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)
1014def _RD(inst, x, y, z, *over):
1015 '''(INTERNAL) Carlson, eqs 2.28 - 2.34.
1016 '''
1017 L = _Lxyz(x, y, z)
1018 a = L.a0(5)
1019 q = L.thresh(_TolRF)
1020 S = _Deferred()
1021 am, m = a, 1
1022 for i in range(_TRIPS):
1023 if fabs(am) > q: # max 7 trips
1024 _iterations(inst, i)
1025 break
1026 t = L[2] # z0...n
1027 a, s, r = L.asr3(a)
1028 S += _1_0 / ((t + r) * s[2] * m)
1029 m *= 4
1030 am = a * m
1031 else: # PYCHOK no cover
1032 raise _no_convergenceError(am, q, Elliptic.fRD, x, y, z, *over,
1033 thresh=True)
1034 x, y = L.rescale(-am, x, y)
1035 xy = x * y
1036 z = (x + y) / _3_0
1037 z2 = z**2
1038 S = S.Fsum.fmul(_3_0)
1039 return _horner(S, am * sqrt(a),
1040 xy - _6_0 * z2,
1041 (xy * _3_0 - _8_0 * z2) * z,
1042 (xy - z2) * _3_0 * z2,
1043 xy * z2 * z, *over)
1046def _RF2(inst, x, y): # 2-arg version, z=0
1047 '''(INTERNAL) Carlson, eqs 2.36 - 2.38.
1048 '''
1049 a, b = sqrt(x), sqrt(y)
1050 if a < b:
1051 a, b = b, a
1052 for i in range(_TRIPS):
1053 t = _TolRG0 * a
1054 if fabs(a - b) <= t: # max 4 trips
1055 _iterations(inst, i)
1056 return (PI / (a + b))
1057 a, b = ((a + b) * _0_5), sqrt(a * b)
1058 else: # PYCHOK no cover
1059 raise _no_convergenceError(a - b, t, Elliptic.fRF, x, y)
1062def _RF3(inst, x, y, z): # 3-arg version
1063 '''(INTERNAL) Carlson, eqs 2.2 - 2.7.
1064 '''
1065 L = _Lxyz(x, y, z)
1066 a = L.a0(3)
1067 q = L.thresh(_TolRF)
1068 am, m = a, 1
1069 for i in range(_TRIPS):
1070 if fabs(am) > q: # max 6 trips
1071 _iterations(inst, i)
1072 break
1073 a, _, _ = L.asr3(a)
1074 m *= 4
1075 am = a * m
1076 else: # PYCHOK no cover
1077 raise _no_convergenceError(am, q, Elliptic.fRF, x, y, z,
1078 thresh=True)
1079 x, y = L.rescale(am, x, y)
1080 z = neg(x + y)
1081 xy = x * y
1082 e2 = xy - z**2
1083 e3 = xy * z
1084 e4 = e2**2
1085 # Polynomial is <https://DLMF.NIST.gov/19.36.E1>
1086 # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44
1087 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16)
1088 # converted to Horner-like form ...
1089 S = Fsum(e4 * 15015, e3 * 6930, e2 * -16380, 17160).fmul(e3)
1090 S += Fsum(e4 * -5775, e2 * 10010, -24024).fmul(e2)
1091 S += Fsum(240240)
1092 return S.fover(sqrt(a) * 240240)
1095def _RG2(inst, x, y): # 2-args and I{doubled}
1096 '''(INTERNAL) Carlson, eqs 2.36 - 2.39.
1097 '''
1098 a, b = sqrt(x), sqrt(y)
1099 if a < b:
1100 a, b = b, a
1101 ab = a - b # fabs(a - b)
1102 S = _Deferred(_0_5 * (a + b)**2)
1103 m = -1
1104 for i in range(_TRIPS): # max 4 trips
1105 t = _TolRG0 * a
1106 if ab <= t:
1107 _iterations(inst, i)
1108 return S.Fsum.fover((a + b) / PI_2)
1109 a, b = ((a + b) * _0_5), sqrt(a * b)
1110 ab = fabs(a - b)
1111 S += ab**2 * m
1112 m *= 2
1113 else: # PYCHOK no cover
1114 raise _no_convergenceError(ab, t, Elliptic.fRG, x, y)
1117def _RG3(inst, x, y, z): # 3-arg version
1118 '''(INTERNAL) Never called with zero B{C{z}}, see C{.fRG}.
1119 '''
1120# if not z:
1121# y, z = z, y
1122 rd = (x - z) * (z - y) # - (y - z)
1123 if rd: # Carlson, eq 1.7
1124 rd = _RD(inst, x, y, z, _3_0 * z / rd)
1125 xyz = x * y
1126 if xyz:
1127 xyz = sqrt(xyz / z**3)
1128 return Fsum(_RF3(inst, x, y, z), rd, xyz).fover(_2_0 / z)
1131def _RJ(inst, x, y, z, p, *over):
1132 '''(INTERNAL) Carlson, eqs 2.17 - 2.25.
1133 '''
1134 def _xyzp(x, y, z, p):
1135 return (x + p) * (y + p) * (z + p)
1137 L = _Lxyz(x, y, z, p)
1138 a = L.a0(5)
1139 q = L.thresh(_TolRD)
1140 S = _Deferred()
1141 n = neg(_xyzp(x, y, z, -p))
1142 am, m = a, 1
1143 for i in range(_TRIPS):
1144 if fabs(am) > q: # max 7 trips
1145 _iterations(inst, i)
1146 break
1147 a, s, _ = L.asr3(a)
1148 d = _xyzp(*s)
1149 if n:
1150 rc = _RC(inst, _1_0, n / d**2 + _1_0)
1151 n *= _1_64th # /= chokes PyChecker
1152 else:
1153 rc = _1_0 # == _RC(None, _1_0, _1_0)
1154 S += rc / (d * m)
1155 m *= 4
1156 am = a * m
1157 else: # PYCHOK no cover
1158 raise _no_convergenceError(am, q, Elliptic.fRJ, x, y, z, p,
1159 thresh=True)
1160 x, y, z = L.rescale(am, x, y, z)
1161 xyz = x * y * z
1162 p = Fsum(x, y, z).fover(_N_2_0)
1163 p2 = p**2
1164 p3 = p2*p
1165 E2 = Fsum(x * y, x * z, y * z, -p2 * _3_0)
1166 E2p = E2 * p
1167 S = S.Fsum.fmul(_6_0)
1168 return _horner(S, am * sqrt(a), E2,
1169 Fsum(p3 * _4_0, xyz, E2p * _2_0),
1170 Fsum(p3 * _3_0, E2p, xyz * _2_0).fmul(p),
1171 xyz * p2, *over)
1173# **) MIT License
1174#
1175# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1176#
1177# Permission is hereby granted, free of charge, to any person obtaining a
1178# copy of this software and associated documentation files (the "Software"),
1179# to deal in the Software without restriction, including without limitation
1180# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1181# and/or sell copies of the Software, and to permit persons to whom the
1182# Software is furnished to do so, subject to the following conditions:
1183#
1184# The above copyright notice and this permission notice shall be included
1185# in all copies or substantial portions of the Software.
1186#
1187# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1188# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1189# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1190# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1191# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1192# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1193# OTHER DEALINGS IN THE SOFTWARE.