Coverage for pygeodesy/elliptic.py: 96%
485 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -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, neg_
78from pygeodesy.constants import EPS, INF, NAN, PI, PI_2, PI_4, \
79 _EPStol as _TolJAC, _0_0, _1_64th, \
80 _0_25, _0_5, _1_0, _2_0, _N_2_0, \
81 _3_0, _4_0, _6_0, _8_0, _180_0, \
82 _360_0, _over
83# from pygeodesy.errors import _ValueError # from .fmath
84from pygeodesy.fmath import fdot, hypot1, zqrt, _ValueError
85from pygeodesy.fsums import Fsum, _sum
86# from pygeodesy.internals import _dunder_nameof # from .lazily
87from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _invalid_, \
88 _invokation_, _negative_, _SPACE_
89from pygeodesy.karney import _K_2_0, _norm180, _signBit, _sincos2
90from pygeodesy.lazily import _ALL_LAZY, _dunder_nameof
91from pygeodesy.named import _Named, _NamedTuple, Fmt, unstr
92from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all
93# from pygeodesy.streprs import Fmt, unstr # from .named
94from pygeodesy.units import Scalar, Scalar_
95# from pygeodesy.utily import sincos2 as _sincos2 # from .karney
97from math import asinh, atan, atan2, ceil, cosh, fabs, floor, \
98 radians, sin, sqrt, tanh
100__all__ = _ALL_LAZY.elliptic
101__version__ = '24.05.18'
103_TolRD = zqrt(EPS * 0.002)
104_TolRF = zqrt(EPS * 0.030)
105_TolRG0 = _TolJAC * 2.7
106_TRIPS = 21 # Max depth, 7 might be sufficient
109class _Cs(object):
110 '''(INTERAL) Complete integrals cache.
111 '''
112 def __init__(self, **kwds):
113 self.__dict__ = kwds
116class _Dsum(list):
117 '''(INTERNAL) Deferred C{Fsum}.
118 '''
119 def __call__(self, s):
120 try: # Fsum *= s
121 return Fsum(*self).fmul(s)
122 except ValueError: # Fsum(NAN) exception
123 return _sum(self) * s
125 def __iadd__(self, x):
126 list.append(self, x)
127 return self
130class Elliptic(_Named):
131 '''Elliptic integrals and functions.
133 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/
134 C++/doc/classGeographicLib_1_1EllipticFunction.html#details>}.
135 '''
136# _alpha2 = 0
137# _alphap2 = 0
138# _eps = EPS
139# _k2 = 0
140# _kp2 = 0
142 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None, **name):
143 '''Constructor, specifying the C{modulus} and C{parameter}.
145 @kwarg name: Optional C{B{name}=NN} (C{str}).
147 @see: Method L{Elliptic.reset} for further details.
149 @note: If only elliptic integrals of the first and second kinds
150 are needed, use C{B{alpha2}=0}, the default value. In
151 that case, we have C{Π(φ, 0, k) = F(φ, k), G(φ, 0, k) =
152 E(φ, k)} and C{H(φ, 0, k) = F(φ, k) - D(φ, k)}.
153 '''
154 self.reset(k2=k2, alpha2=alpha2, kp2=kp2, alphap2=alphap2)
155 if name:
156 self.name = name
158 @Property_RO
159 def alpha2(self):
160 '''Get α^2, the square of the parameter (C{float}).
161 '''
162 return self._alpha2
164 @Property_RO
165 def alphap2(self):
166 '''Get α'^2, the square of the complementary parameter (C{float}).
167 '''
168 return self._alphap2
170 @Property_RO
171 def cD(self):
172 '''Get Jahnke's complete integral C{D(k)} (C{float}),
173 U{defined<https://DLMF.NIST.gov/19.2.E6>}.
174 '''
175 return self._cDEKEeps.cD
177 @Property_RO
178 def _cDEKEeps(self):
179 '''(INTERNAL) Get the complete integrals D, E, K and KE plus C{eps}.
180 '''
181 k2, kp2 = self.k2, self.kp2
182 if k2:
183 if kp2:
184 try:
185 self._iteration = 0
186 # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3
187 # <https://DLMF.NIST.gov/19.25.E1>
188 D = _RD(self, _0_0, kp2, _1_0, _3_0)
189 cD = float(D)
190 # Complete elliptic integral E(k), Carlson eq. 4.2
191 # <https://DLMF.NIST.gov/19.25.E1>
192 cE = _rG2(self, kp2, _1_0, PI_=PI_2)
193 # Complete elliptic integral K(k), Carlson eq. 4.1
194 # <https://DLMF.NIST.gov/19.25.E1>
195 cK = _rF2(self, kp2, _1_0)
196 cKE = float(D.fmul(k2))
197 eps = k2 / (sqrt(kp2) + _1_0)**2
199 except Exception as e:
200 raise _ellipticError(self.reset, k2=k2, kp2=kp2, cause=e)
201 else:
202 cD = cK = cKE = INF
203 cE = _1_0
204 eps = k2
205 else:
206 cD = PI_4
207 cE = cK = PI_2
208 cKE = _0_0 # k2 * cD
209 eps = EPS
211 return _Cs(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps)
213 @Property_RO
214 def cE(self):
215 '''Get the complete integral of the second kind C{E(k)}
216 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}.
217 '''
218 return self._cDEKEeps.cE
220 @Property_RO
221 def cG(self):
222 '''Get Legendre's complete geodesic longitude integral
223 C{G(α^2, k)} (C{float}).
224 '''
225 return self._cGHPi.cG
227 @Property_RO
228 def _cGHPi(self):
229 '''(INTERNAL) Get the complete integrals G, H and Pi.
230 '''
231 alpha2, alphap2, kp2 = self.alpha2, self.alphap2, self.kp2
232 try:
233 self._iteration = 0
234 if alpha2:
235 if alphap2:
236 if kp2: # <https://DLMF.NIST.gov/19.25.E2>
237 cK = self.cK
238 Rj = _RJ(self, _0_0, kp2, _1_0, alphap2, _3_0)
239 cG = float(Rj * (alpha2 - self.k2) + cK) # G(alpha2, k)
240 cH = -float(Rj * alphap2 - cK) # H(alpha2, k)
241 cPi = float(Rj * alpha2 + cK) # Pi(alpha2, k)
242 else:
243 cG = cH = _rC(self, _1_0, alphap2)
244 cPi = INF # XXX or NAN?
245 else:
246 cG = cH = cPi = INF # XXX or NAN?
247 else:
248 cG, cPi = self.cE, self.cK
249 # H = K - D but this involves large cancellations if k2 is near 1.
250 # So write (for alpha2 = 0)
251 # H = int(cos(phi)**2 / sqrt(1-k2 * sin(phi)**2), phi, 0, pi/2)
252 # = 1 / sqrt(1-k2) * int(sin(phi)**2 / sqrt(1-k2/kp2 * sin(phi)**2,...)
253 # = 1 / kp * D(i * k/kp)
254 # and use D(k) = RD(0, kp2, 1) / 3, so
255 # H = 1/kp * RD(0, 1/kp2, 1) / 3
256 # = kp2 * RD(0, 1, kp2) / 3
257 # using <https://DLMF.NIST.gov/19.20.E18>. Equivalently
258 # RF(x, 1) - RD(0, x, 1) / 3 = x * RD(0, 1, x) / 3 for x > 0
259 # For k2 = 1 and alpha2 = 0, we have
260 # H = int(cos(phi),...) = 1
261 cH = float(_RD(self, _0_0, _1_0, kp2, _3_0 / kp2)) if kp2 else _1_0
263 except Exception as e:
264 raise _ellipticError(self.reset, kp2=kp2, alpha2 =alpha2,
265 alphap2=alphap2, cause=e)
266 return _Cs(cG=cG, cH=cH, cPi=cPi)
268 @Property_RO
269 def cH(self):
270 '''Get Cayley's complete geodesic longitude difference integral
271 C{H(α^2, k)} (C{float}).
272 '''
273 return self._cGHPi.cH
275 @Property_RO
276 def cK(self):
277 '''Get the complete integral of the first kind C{K(k)}
278 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}.
279 '''
280 return self._cDEKEeps.cK
282 @Property_RO
283 def cKE(self):
284 '''Get the difference between the complete integrals of the
285 first and second kinds, C{K(k) − E(k)} (C{float}).
286 '''
287 return self._cDEKEeps.cKE
289 @Property_RO
290 def cPi(self):
291 '''Get the complete integral of the third kind C{Pi(α^2, k)}
292 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}.
293 '''
294 return self._cGHPi.cPi
296 def deltaD(self, sn, cn, dn):
297 '''Jahnke's periodic incomplete elliptic integral.
299 @arg sn: sin(φ).
300 @arg cn: cos(φ).
301 @arg dn: sqrt(1 − k2 * sin(2φ)).
303 @return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}).
305 @raise EllipticError: Invalid invokation or no convergence.
306 '''
307 return _deltaX(sn, cn, dn, self.cD, self.fD)
309 def deltaE(self, sn, cn, dn):
310 '''The periodic incomplete integral of the second kind.
312 @arg sn: sin(φ).
313 @arg cn: cos(φ).
314 @arg dn: sqrt(1 − k2 * sin(2φ)).
316 @return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}).
318 @raise EllipticError: Invalid invokation or no convergence.
319 '''
320 return _deltaX(sn, cn, dn, self.cE, self.fE)
322 def deltaEinv(self, stau, ctau):
323 '''The periodic inverse of the incomplete integral of the second kind.
325 @arg stau: sin(τ)
326 @arg ctau: cos(τ)
328 @return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}).
330 @raise EllipticError: No convergence.
331 '''
332 try:
333 if _signBit(ctau): # pi periodic
334 stau, ctau = neg_(stau, ctau)
335 t = atan2(stau, ctau)
336 return self._Einv(t * self.cE / PI_2) - t
338 except Exception as e:
339 raise _ellipticError(self.deltaEinv, stau, ctau, cause=e)
341 def deltaF(self, sn, cn, dn):
342 '''The periodic incomplete integral of the first kind.
344 @arg sn: sin(φ).
345 @arg cn: cos(φ).
346 @arg dn: sqrt(1 − k2 * sin(2φ)).
348 @return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}).
350 @raise EllipticError: Invalid invokation or no convergence.
351 '''
352 return _deltaX(sn, cn, dn, self.cK, self.fF)
354 def deltaG(self, sn, cn, dn):
355 '''Legendre's periodic geodesic longitude integral.
357 @arg sn: sin(φ).
358 @arg cn: cos(φ).
359 @arg dn: sqrt(1 − k2 * sin(2φ)).
361 @return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}).
363 @raise EllipticError: Invalid invokation or no convergence.
364 '''
365 return _deltaX(sn, cn, dn, self.cG, self.fG)
367 def deltaH(self, sn, cn, dn):
368 '''Cayley's periodic geodesic longitude difference integral.
370 @arg sn: sin(φ).
371 @arg cn: cos(φ).
372 @arg dn: sqrt(1 − k2 * sin(2φ)).
374 @return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}).
376 @raise EllipticError: Invalid invokation or no convergence.
377 '''
378 return _deltaX(sn, cn, dn, self.cH, self.fH)
380 def deltaPi(self, sn, cn, dn):
381 '''The periodic incomplete integral of the third kind.
383 @arg sn: sin(φ).
384 @arg cn: cos(φ).
385 @arg dn: sqrt(1 − k2 * sin(2φ)).
387 @return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ
388 (C{float}).
390 @raise EllipticError: Invalid invokation or no convergence.
391 '''
392 return _deltaX(sn, cn, dn, self.cPi, self.fPi)
394 def _Einv(self, x):
395 '''(INTERNAL) Helper for C{.deltaEinv} and C{.fEinv}.
396 '''
397 E2 = self.cE * _2_0
398 n = floor(x / E2 + _0_5)
399 r = x - E2 * n # r in [-cE, cE)
400 # linear approximation
401 phi = PI * r / E2 # phi in [-PI_2, PI_2)
402 Phi = Fsum(phi)
403 # first order correction
404 phi = Phi.fsum_(self.eps * sin(phi * _2_0) / _N_2_0)
405 self._iteration = 0
406 # For kp2 close to zero use asin(r / cE) or J. P. Boyd,
407 # Applied Math. and Computation 218, 7005-7013 (2012)
408 # <https://DOI.org/10.1016/j.amc.2011.12.021>
409 _Phi2 = Phi.fsum2f_ # aggregate
410 _abs = fabs
411 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC
412 sn, cn, dn = self._sncndn3(phi)
413 if dn:
414 sn = self.fE(sn, cn, dn)
415 phi, d = _Phi2((r - sn) / dn)
416 else: # PYCHOK no cover
417 d = _0_0 # XXX continue?
418 if _abs(d) < _TolJAC: # 3-4 trips
419 _iterations(self, i)
420 break
421 else: # PYCHOK no cover
422 raise _convergenceError(d, _TolJAC)
423 return Phi.fsum_(n * PI) if n else phi
425 @Property_RO
426 def eps(self):
427 '''Get epsilon (C{float}).
428 '''
429 return self._cDEKEeps.eps
431 def fD(self, phi_or_sn, cn=None, dn=None):
432 '''Jahnke's incomplete elliptic integral in terms of
433 Jacobi elliptic functions.
435 @arg phi_or_sn: φ or sin(φ).
436 @kwarg cn: C{None} or cos(φ).
437 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
439 @return: D(φ, k) as though φ ∈ (−π, π] (C{float}),
440 U{defined<https://DLMF.NIST.gov/19.2.E4>}.
442 @raise EllipticError: Invalid invokation or no convergence.
443 '''
444 def _fD(sn, cn, dn):
445 r = fabs(sn)**3
446 if r:
447 r = float(_RD(self, cn**2, dn**2, _1_0, _3_0 / r))
448 return r
450 return self._fXf(phi_or_sn, cn, dn, self.cD,
451 self.deltaD, _fD)
453 def fDelta(self, sn, cn):
454 '''The C{Delta} amplitude function.
456 @arg sn: sin(φ).
457 @arg cn: cos(φ).
459 @return: sqrt(1 − k2 * sin(2φ)) (C{float}).
460 '''
461 try:
462 k2 = self.k2
463 s = (self.kp2 + cn**2 * k2) if k2 > 0 else (
464 (_1_0 - sn**2 * k2) if k2 < 0 else self.kp2)
465 return sqrt(s) if s else _0_0
467 except Exception as e:
468 raise _ellipticError(self.fDelta, sn, cn, k2=k2, cause=e)
470 def fE(self, phi_or_sn, cn=None, dn=None):
471 '''The incomplete integral of the second kind in terms of
472 Jacobi elliptic functions.
474 @arg phi_or_sn: φ or sin(φ).
475 @kwarg cn: C{None} or cos(φ).
476 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
478 @return: E(φ, k) as though φ ∈ (−π, π] (C{float}),
479 U{defined<https://DLMF.NIST.gov/19.2.E5>}.
481 @raise EllipticError: Invalid invokation or no convergence.
482 '''
483 def _fE(sn, cn, dn):
484 '''(INTERNAL) Core of C{.fE}.
485 '''
486 if sn:
487 sn2, cn2, dn2 = sn**2, cn**2, dn**2
488 kp2, k2 = self.kp2, self.k2
489 if k2 <= 0: # Carlson, eq. 4.6, <https://DLMF.NIST.gov/19.25.E9>
490 Ei = _RF3(self, cn2, dn2, _1_0)
491 if k2:
492 Ei -= _RD(self, cn2, dn2, _1_0, _3over(k2, sn2))
493 elif kp2 >= 0: # k2 > 0, <https://DLMF.NIST.gov/19.25.E10>
494 Ei = _over(k2 * fabs(cn), dn) # float
495 if kp2:
496 Ei += (_RD( self, cn2, _1_0, dn2, _3over(k2, sn2)) +
497 _RF3(self, cn2, dn2, _1_0)) * kp2
498 else: # kp2 < 0, <https://DLMF.NIST.gov/19.25.E11>
499 Ei = _over(dn, fabs(cn))
500 Ei -= _RD(self, dn2, _1_0, cn2, _3over(kp2, sn2))
501 Ei *= fabs(sn)
502 ei = float(Ei)
503 else: # PYCHOK no cover
504 ei = _0_0
505 return ei
507 return self._fXf(phi_or_sn, cn, dn, self.cE,
508 self.deltaE, _fE)
510 def fEd(self, deg):
511 '''The incomplete integral of the second kind with
512 the argument given in C{degrees}.
514 @arg deg: Angle (C{degrees}).
516 @return: E(π B{C{deg}} / 180, k) (C{float}).
518 @raise EllipticError: No convergence.
519 '''
520 if _K_2_0:
521 e = round((deg - _norm180(deg)) / _360_0)
522 elif fabs(deg) < _180_0:
523 e = _0_0
524 else:
525 e = ceil(deg / _360_0 - _0_5)
526 deg -= e * _360_0
527 return self.fE(radians(deg)) + e * self.cE * _4_0
529 def fEinv(self, x):
530 '''The inverse of the incomplete integral of the second kind.
532 @arg x: Argument (C{float}).
534 @return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}}
535 (C{float}).
537 @raise EllipticError: No convergence.
538 '''
539 try:
540 return self._Einv(x)
541 except Exception as e:
542 raise _ellipticError(self.fEinv, x, cause=e)
544 def fF(self, phi_or_sn, cn=None, dn=None):
545 '''The incomplete integral of the first kind in terms of
546 Jacobi elliptic functions.
548 @arg phi_or_sn: φ or sin(φ).
549 @kwarg cn: C{None} or cos(φ).
550 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
552 @return: F(φ, k) as though φ ∈ (−π, π] (C{float}),
553 U{defined<https://DLMF.NIST.gov/19.2.E4>}.
555 @raise EllipticError: Invalid invokation or no convergence.
556 '''
557 def _fF(sn, cn, dn):
558 r = fabs(sn)
559 if r:
560 r = float(_RF3(self, cn**2, dn**2, _1_0).fmul(r))
561 return r
563 return self._fXf(phi_or_sn, cn, dn, self.cK,
564 self.deltaF, _fF)
566 def fG(self, phi_or_sn, cn=None, dn=None):
567 '''Legendre's geodesic longitude integral 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: G(φ, k) as though φ ∈ (−π, π] (C{float}).
576 @raise EllipticError: Invalid invokation or no convergence.
578 @note: Legendre expresses the longitude of a point on the
579 geodesic in terms of this combination of elliptic
580 integrals in U{Exercices de Calcul Intégral, Vol 1
581 (1811), p 181<https://Books.Google.com/books?id=
582 riIOAAAAQAAJ&pg=PA181>}.
584 @see: U{Geodesics in terms of elliptic integrals<https://
585 GeographicLib.SourceForge.io/C++/doc/geodesic.html#geodellip>}
586 for the expression for the longitude in terms of this function.
587 '''
588 return self._fXa(phi_or_sn, cn, dn, self.alpha2 - self.k2,
589 self.cG, self.deltaG)
591 def fH(self, phi_or_sn, cn=None, dn=None):
592 '''Cayley's geodesic longitude difference integral in terms of
593 Jacobi elliptic functions.
595 @arg phi_or_sn: φ or sin(φ).
596 @kwarg cn: C{None} or cos(φ).
597 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
599 @return: H(φ, k) as though φ ∈ (−π, π] (C{float}).
601 @raise EllipticError: Invalid invokation or no convergence.
603 @note: Cayley expresses the longitude difference of a point
604 on the geodesic in terms of this combination of
605 elliptic integrals in U{Phil. Mag. B{40} (1870), p 333
606 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}.
608 @see: U{Geodesics in terms of elliptic integrals<https://
609 GeographicLib.SourceForge.io/C++/doc/geodesic.html#geodellip>}
610 for the expression for the longitude in terms of this function.
611 '''
612 return self._fXa(phi_or_sn, cn, dn, -self.alphap2,
613 self.cH, self.deltaH)
615 def fPi(self, phi_or_sn, cn=None, dn=None):
616 '''The incomplete integral of the third kind in terms of
617 Jacobi elliptic functions.
619 @arg phi_or_sn: φ or sin(φ).
620 @kwarg cn: C{None} or cos(φ).
621 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)).
623 @return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}).
625 @raise EllipticError: Invalid invokation or no convergence.
626 '''
627 if dn is None and cn is not None: # and isscalar(phi_or_sn)
628 dn = self.fDelta(phi_or_sn, cn) # in .triaxial
629 return self._fXa(phi_or_sn, cn, dn, self.alpha2,
630 self.cPi, self.deltaPi)
632 def _fXa(self, phi_or_sn, cn, dn, aX, cX, deltaX):
633 '''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}.
634 '''
635 def _fX(sn, cn, dn):
636 if sn:
637 cn2, dn2 = cn**2, dn**2
638 R = _RF3(self, cn2, dn2, _1_0)
639 if aX:
640 sn2 = sn**2
641 p = sn2 * self.alphap2 + cn2
642 R += _RJ(self, cn2, dn2, _1_0, p, _3over(aX, sn2))
643 R *= fabs(sn)
644 r = float(R)
645 else: # PYCHOK no cover
646 r = _0_0
647 return r
649 return self._fXf(phi_or_sn, cn, dn, cX, deltaX, _fX)
651 def _fXf(self, phi_or_sn, cn, dn, cX, deltaX, fX):
652 '''(INTERNAL) Helper for C{.fD}, C{.fE}, C{.fF} and C{._fXa}.
653 '''
654 self._iteration = 0 # aggregate
655 phi = sn = phi_or_sn
656 if cn is dn is None: # fX(phi) call
657 sn, cn, dn = self._sncndn3(phi)
658 if fabs(phi) >= PI:
659 return (deltaX(sn, cn, dn) + phi) * cX / PI_2
660 # fall through
661 elif cn is None or dn is None:
662 n = NN(_f_, deltaX.__name__[5:])
663 raise _ellipticError(n, sn, cn, dn)
665 if _signBit(cn): # enforce usual trig-like symmetries
666 xi = cX * _2_0 - fX(sn, cn, dn)
667 else:
668 xi = fX(sn, cn, dn) if cn > 0 else cX
669 return copysign0(xi, sn)
671 @Property_RO
672 def k2(self):
673 '''Get k^2, the square of the modulus (C{float}).
674 '''
675 return self._k2
677 @Property_RO
678 def kp2(self):
679 '''Get k'^2, the square of the complementary modulus (C{float}).
680 '''
681 return self._kp2
683 def reset(self, k2=0, alpha2=0, kp2=None, alphap2=None): # MCCABE 13
684 '''Reset the modulus, parameter and the complementaries.
686 @kwarg k2: Modulus squared (C{float}, NINF <= k^2 <= 1).
687 @kwarg alpha2: Parameter squared (C{float}, NINF <= α^2 <= 1).
688 @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0).
689 @kwarg alphap2: Complementary parameter squared (C{float}, α'^2 >= 0).
691 @raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}}
692 or B{C{alphap2}}.
694 @note: The arguments must satisfy C{B{k2} + B{kp2} = 1} and
695 C{B{alpha2} + B{alphap2} = 1}. No checking is done
696 that these conditions are met to enable accuracy to be
697 maintained, e.g., when C{k} is very close to unity.
698 '''
699 if self.__dict__:
700 _update_all(self, _Named.iteration._uname, Base=Property_RO)
702 self._k2 = Scalar_(k2=k2, Error=EllipticError, low=None, high=_1_0)
703 self._kp2 = Scalar_(kp2=((_1_0 - k2) if kp2 is None else kp2), Error=EllipticError)
705 self._alpha2 = Scalar_(alpha2=alpha2, Error=EllipticError, low=None, high=_1_0)
706 self._alphap2 = Scalar_(alphap2=((_1_0 - alpha2) if alphap2 is None else alphap2),
707 Error=EllipticError)
709 # Values of complete elliptic integrals for k = 0,1 and alpha = 0,1
710 # K E D
711 # k = 0: pi/2 pi/2 pi/4
712 # k = 1: inf 1 inf
713 # Pi G H
714 # k = 0, alpha = 0: pi/2 pi/2 pi/4
715 # k = 1, alpha = 0: inf 1 1
716 # k = 0, alpha = 1: inf inf pi/2
717 # k = 1, alpha = 1: inf inf inf
718 #
719 # G(0, k) = Pi(0, k) = H(1, k) = E(k)
720 # H(0, k) = K(k) - D(k)
721 # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2))
722 # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1))
723 # Pi(alpha2, 1) = inf
724 # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2)
726 def sncndn(self, x):
727 '''The Jacobi elliptic function.
729 @arg x: The argument (C{float}).
731 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with
732 C{*n(B{x}, k)}.
734 @raise EllipticError: No convergence.
735 '''
736 self._iteration = 0 # reset
737 try: # Bulirsch's sncndn routine, p 89.
738 if self.kp2:
739 c, d, cd, mn = self._sncndn4
740 dn = _1_0
741 sn, cn = _sincos2(x * cd)
742 if sn:
743 a = cn / sn
744 c *= a
745 for m, n in reversed(mn):
746 a *= c
747 c *= dn
748 dn = (n + a) / (m + a)
749 a = c / m
750 a = _1_0 / hypot1(c)
751 sn = neg(a) if _signBit(sn) else a
752 cn = c * sn
753 if d and _signBit(self.kp2):
754 cn, dn = dn, cn
755 sn = sn / d # /= chokes PyChecker
756 else:
757 sn = tanh(x)
758 cn = dn = _1_0 / cosh(x)
760 except Exception as e:
761 raise _ellipticError(self.sncndn, x, kp2=self.kp2, cause=e)
763 return Elliptic3Tuple(sn, cn, dn, iteration=self._iteration)
765 def _sncndn3(self, phi):
766 '''(INTERNAL) Helper for C{.fEinv} and C{._fXf}.
767 '''
768 sn, cn = _sincos2(phi)
769 return sn, cn, self.fDelta(sn, cn)
771 @Property_RO
772 def _sncndn4(self):
773 '''(INTERNAL) Get Bulirsch' 4-tuple C{(c, d, cd, mn)}.
774 '''
775 # Bulirsch's sncndn routine, p 89.
776 d, mc = 0, self.kp2
777 if _signBit(mc):
778 d = _1_0 - mc
779 mc = neg(mc / d)
780 d = sqrt(d)
782 mn, a = [], _1_0
783 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC
784 mc = sqrt(mc)
785 mn.append((a, mc))
786 c = (a + mc) * _0_5
787 r = fabs(mc - a)
788 t = _TolJAC * a
789 if r <= t: # 6 trips, quadratic
790 _iterations(self, i)
791 break
792 mc *= a
793 a = c
794 else: # PYCHOK no cover
795 raise _convergenceError(r, t)
796 cd = (c * d) if d else c
797 return c, d, cd, mn
799 @staticmethod
800 def fRC(x, y):
801 '''Degenerate symmetric integral of the first kind C{RC(x, y)}.
803 @return: C{RC(x, y)}, equivalent to C{RF(x, y, y)}.
805 @see: U{C{RC} definition<https://DLMF.NIST.gov/19.2.E17>} and
806 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
807 '''
808 return _rC(None, x, y)
810 @staticmethod
811 def fRD(x, y, z, *over):
812 '''Degenerate symmetric integral of the third kind C{RD(x, y, z)}.
814 @return: C{RD(x, y, z) / over}, equivalent to C{RJ(x, y, z, z)
815 / over} with C{over} typically 3.
817 @see: U{C{RD} definition<https://DLMF.NIST.gov/19.16.E5>} and
818 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
819 '''
820 try:
821 return float(_RD(None, x, y, z, *over))
822 except Exception as e:
823 raise _ellipticError(Elliptic.fRD, x, y, z, *over, cause=e)
825 @staticmethod
826 def fRF(x, y, z=0):
827 '''Symmetric or complete symmetric integral of the first kind
828 C{RF(x, y, z)} respectively C{RF(x, y)}.
830 @return: C{RF(x, y, z)} or C{RF(x, y)} for missing or zero B{C{z}}.
832 @see: U{C{RF} definition<https://DLMF.NIST.gov/19.16.E1>} and
833 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
834 '''
835 try:
836 return float(_RF3(None, x, y, z)) if z else _rF2(None, x, y)
837 except Exception as e:
838 raise _ellipticError(Elliptic.fRF, x, y, z, cause=e)
840 @staticmethod
841 def fRG(x, y, z=0):
842 '''Symmetric or complete symmetric integral of the second kind
843 C{RG(x, y, z)} respectively C{RG(x, y)}.
845 @return: C{RG(x, y, z)} or C{RG(x, y)} for missing or zero B{C{z}}.
847 @see: U{C{RG} definition<https://DLMF.NIST.gov/19.16.E3>},
848 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>} and
849 U{RG<https://GeographicLib.SourceForge.io/C++/doc/
850 EllipticFunction_8cpp_source.html#l00096>} version 2.3.
851 '''
852 try:
853 return _rG2(None, x, y) if z == 0 else (
854 _rG2(None, z, x) if y == 0 else (
855 _rG2(None, y, z) if x == 0 else _rG3(None, x, y, z)))
856 except Exception as e:
857 t = _negative_ if min(x, y, z) < 0 else NN
858 raise _ellipticError(Elliptic.fRG, x, y, z, cause=e, txt=t)
860 @staticmethod
861 def fRJ(x, y, z, p):
862 '''Symmetric integral of the third kind C{RJ(x, y, z, p)}.
864 @return: C{RJ(x, y, z, p)}.
866 @see: U{C{RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and
867 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.
868 '''
869 try:
870 return float(_RJ(None, x, y, z, p))
871 except Exception as e:
872 raise _ellipticError(Elliptic.fRJ, x, y, z, p, cause=e)
874 @staticmethod
875 def _RFRD(x, y, z, m):
876 # in .auxilats.AuxDLat.DE, .auxilats.AuxLat.Rectifying
877 try: # float(RF(x, y, z) - RD(x, y, z, 3 / m))
878 R = _RF3(None, x, y, z)
879 if m:
880 R -= _RD(None, x, y, z, _3_0 / m)
881 except Exception as e:
882 raise _ellipticError(Elliptic._RFRD, x, y, z, m, cause=e)
883 return float(R)
885_allPropertiesOf_n(15, Elliptic) # # PYCHOK assert, see Elliptic.reset
888class EllipticError(_ValueError):
889 '''Elliptic function, integral, convergence or other L{Elliptic} issue.
890 '''
891 pass
894class Elliptic3Tuple(_NamedTuple):
895 '''3-Tuple C{(sn, cn, dn)} all C{scalar}.
896 '''
897 _Names_ = ('sn', 'cn', 'dn')
898 _Units_ = ( Scalar, Scalar, Scalar)
901class _List(list):
902 '''(INTERNAL) Helper for C{_RD}, C{_RF3} and C{_RJ}.
903 '''
904 _a0 = None
905# _xyzp = ()
907 def __init__(self, *xyzp): # x, y, z [, p]
908 list.__init__(self, xyzp)
909 self._xyzp = xyzp
911 def a0(self, n):
912 '''Compute the initial C{a}.
913 '''
914 t = tuple(self)
915 m = n - len(t)
916 if m > 0:
917 t += t[-1:] * m
918 try:
919 a = Fsum(*t).fover(n)
920 except ValueError: # Fsum(NAN) exception
921 a = _sum(t) / n
922 self._a0 = a
923 return a
925 def amrs4(self, inst, y, Tol):
926 '''Yield Carlson 4-tuples C{(An, mul, lam, s)} plus sentinel, with
927 C{lam = fdot(sqrt(x), ... (z))} and C{s = (sqrt(x), ... (p))}.
928 '''
929 L = self
930 a = L.a0(5 if y else 3)
931 m = 1
932 t = max(fabs(a - _) for _ in L) / Tol
933 for i in range(_TRIPS):
934 d = fabs(a * m)
935 if d > t: # 3-6 trips
936 _iterations(inst, i)
937 break
938 s = map2(sqrt, L) # sqrt(x), srqt(y), sqrt(z) [, sqrt(p)]
939 try:
940 r = fdot(s[:3], s[1], s[2], s[0]) # sqrt(x) * sqrt(y) + ...
941 except ValueError: # Fsum(NAN) exception
942 r = _sum(s[i] * s[(i + 1) % 3] for i in range(3))
943 L[:] = ((r + _) * _0_25 for _ in L)
944 a = (r + a) * _0_25
945 if y: # yield only if used
946 yield a, m, r, s # L[2] is next z
947 m *= 4
948 else: # PYCHOK no cover
949 raise _convergenceError(d, t, thresh=True)
950 yield a, m, None, () # sentinel: same a, next m, no r and s
952 def rescale(self, am, *xs):
953 '''Rescale C{x}, C{y}, ...
954 '''
955 # assert am
956 a0 = self._a0
957 _am = _1_0 / am
958 for x in xs:
959 yield (a0 - x) * _am
962def _ab2(inst, x, y):
963 '''(INTERNAL) Yield Carlson 2-tuples C{(xn, yn)}.
964 '''
965 a, b = sqrt(x), sqrt(y)
966 if b > a:
967 b, a = a, b
968 for i in range(_TRIPS):
969 yield a, b # xi, yi
970 d = fabs(a - b)
971 t = _TolRG0 * a
972 if d <= t: # 3-4 trips
973 _iterations(inst, i)
974 break
975 a, b = ((a + b) * _0_5), sqrt(a * b)
976 else: # PYCHOK no cover
977 raise _convergenceError(d, t)
980def _convergenceError(d, tol, **thresh):
981 '''(INTERNAL) Format a no-convergence Error.
982 '''
983 t = Fmt.no_convergence(d, tol, **thresh)
984 return ValueError(t) # txt only
987def _deltaX(sn, cn, dn, cX, fX):
988 '''(INTERNAL) Helper for C{Elliptic.deltaD} thru C{.deltaPi}.
989 '''
990 try:
991 if cn is None or dn is None:
992 raise ValueError(_invalid_)
994 if _signBit(cn):
995 sn, cn = neg_(sn, cn)
996 r = fX(sn, cn, dn) * PI_2 / cX
997 return r - atan2(sn, cn)
999 except Exception as e:
1000 n = NN(_delta_, fX.__name__[1:])
1001 raise _ellipticError(n, sn, cn, dn, cause=e)
1004def _ellipticError(where, *args, **kwds_cause_txt):
1005 '''(INTERNAL) Format an L{EllipticError}.
1006 '''
1007 def _x_t_kwds(cause=None, txt=NN, **kwds):
1008 return cause, txt, kwds
1010 x, t, kwds = _x_t_kwds(**kwds_cause_txt)
1012 n = _dunder_nameof(where, where)
1013 n = _DOT_(Elliptic.__name__, n)
1014 n = _SPACE_(_invokation_, n)
1015 u = unstr(n, *args, **kwds)
1016 return EllipticError(u, cause=x, txt=t)
1019def _Horner(S, e1, E2, E3, E4, E5, *over):
1020 '''(INTERNAL) Horner form for C{_RD} and C{_RJ} below.
1021 '''
1022 E22 = E2**2
1023 # Polynomial is <https://DLMF.NIST.gov/19.36.E2>
1024 # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52
1025 # + 3*E5/26 - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20
1026 # + 45*E2**2*E3/272 - 9*(E3*E4+E2*E5)/68)
1027 # converted to Horner-like form ...
1028 e = e1 * 4084080
1029 S *= e
1030 S += Fsum(E2 * -540540, 471240).fmul(E5)
1031 S += Fsum(E2 * 612612, E3 * -540540, -556920).fmul(E4)
1032 S += Fsum(E2 * -706860, E22 * 675675, E3 * 306306, 680680).fmul(E3)
1033 S += Fsum(E2 * 417690, E22 * -255255, -875160).fmul(E2)
1034 S += 4084080
1035 if over:
1036 e *= over[0]
1037 return S.fdiv(e) # Fsum
1040def _iterations(inst, i):
1041 '''(INTERNAL) Aggregate iterations B{C{i}}.
1042 '''
1043 if inst and i > 0:
1044 inst._iteration += i
1047def _3over(a, b):
1048 '''(INTERNAL) Return C{3 / (a * b)}.
1049 '''
1050 return _over(_3_0, a * b)
1053def _rC(unused, x, y):
1054 '''(INTERNAL) Defined only for C{y != 0} and C{x >= 0}.
1055 '''
1056 d = x - y
1057 if d < 0: # catch NaN
1058 # <https://DLMF.NIST.gov/19.2.E18>
1059 d = -d
1060 r = atan(sqrt(d / x)) if x > 0 else PI_2
1061 elif d == 0: # XXX d < EPS0? or EPS02 or _EPSmin
1062 d, r = y, _1_0
1063 elif y > 0: # <https://DLMF.NIST.gov/19.2.E19>
1064 r = asinh(sqrt(d / y)) # atanh(sqrt((x - y) / x))
1065 elif y < 0: # <https://DLMF.NIST.gov/19.2.E20>
1066 r = asinh(sqrt(-x / y)) # atanh(sqrt(x / (x - y)))
1067 else: # PYCHOK no cover
1068 raise _ellipticError(Elliptic.fRC, x, y)
1069 return r / sqrt(d) # float
1072def _RD(inst, x, y, z, *over):
1073 '''(INTERNAL) Carlson, eqs 2.28 - 2.34.
1074 '''
1075 L = _List(x, y, z)
1076 S = _Dsum()
1077 for a, m, r, s in L.amrs4(inst, True, _TolRF):
1078 if s:
1079 S += _over(_3_0, (r + z) * s[2] * m)
1080 z = L[2] # s[2] = sqrt(z)
1081 x, y = L.rescale(-a * m, x, y)
1082 xy = x * y
1083 z = (x + y) / _3_0
1084 z2 = z**2
1085 return _Horner(S(_1_0), sqrt(a) * a * m,
1086 xy - _6_0 * z2,
1087 (xy * _3_0 - _8_0 * z2) * z,
1088 (xy - z2) * _3_0 * z2,
1089 xy * z2 * z, *over) # Fsum
1092def _rF2(inst, x, y): # 2-arg version, z=0
1093 '''(INTERNAL) Carlson, eqs 2.36 - 2.38.
1094 '''
1095 for a, b in _ab2(inst, x, y): # PYCHOK yield
1096 pass
1097 return _over(PI, a + b) # float
1100def _RF3(inst, x, y, z): # 3-arg version
1101 '''(INTERNAL) Carlson, eqs 2.2 - 2.7.
1102 '''
1103 L = _List(x, y, z)
1104 for a, m, _, _ in L.amrs4(inst, False, _TolRF):
1105 pass
1106 x, y = L.rescale(a * m, x, y)
1107 z = neg(x + y)
1108 xy = x * y
1109 e2 = xy - z**2
1110 e3 = xy * z
1111 e4 = e2**2
1112 # Polynomial is <https://DLMF.NIST.gov/19.36.E1>
1113 # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44
1114 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16)
1115 # converted to Horner-like form ...
1116 S = Fsum(e4 * 15015, e3 * 6930, e2 * -16380, 17160).fmul(e3)
1117 S += Fsum(e4 * -5775, e2 * 10010, -24024).fmul(e2)
1118 S += 240240
1119 return S.fdiv(sqrt(a) * 240240) # Fsum
1122def _rG2(inst, x, y, PI_=PI_4): # 2-args
1123 '''(INTERNAL) Carlson, eqs 2.36 - 2.39.
1124 '''
1125 m = -1 # neg!
1126 S = None
1127 for a, b in _ab2(inst, x, y): # PYCHOK yield
1128 if S is None: # initial
1129 S = _Dsum()
1130 S += (a + b)**2 * _0_5
1131 else:
1132 S += (a - b)**2 * m
1133 m *= 2
1134 return S(PI_).fover(a + b)
1137def _rG3(inst, x, y, z): # 3-arg version
1138 '''(INTERNAL) C{x}, C{y} and C{z} all non-zero, see C{.fRG}.
1139 '''
1140 R = _RF3(inst, x, y, z) * z
1141 rd = (x - z) * (z - y) # - (y - z)
1142 if rd: # Carlson, eq 1.7
1143 R += _RD(inst, x, y, z, _3_0 / rd)
1144 R += sqrt(x * y / z)
1145 return R.fover(_2_0)
1148def _RJ(inst, x, y, z, p, *over):
1149 '''(INTERNAL) Carlson, eqs 2.17 - 2.25.
1150 '''
1151 def _xyzp(x, y, z, p):
1152 return (x + p) * (y + p) * (z + p)
1154 L = _List(x, y, z, p)
1155 n = neg(_xyzp(x, y, z, -p))
1156 S = _Dsum()
1157 for a, m, _, s in L.amrs4(inst, True, _TolRD):
1158 if s:
1159 d = _xyzp(*s)
1160 if d:
1161 if n:
1162 rc = _rC(inst, _1_0, n / d**2 + _1_0)
1163 n *= _1_64th # /= chokes PyChecker
1164 else:
1165 rc = _1_0 # == _rC(None, _1_0, _1_0)
1166 S += rc / (d * m)
1167 else: # PYCHOK no cover
1168 return NAN
1169 x, y, z = L.rescale(a * m, x, y, z)
1170 p = Fsum(x, y, z).fover(_N_2_0)
1171 p2 = p**2
1172 p3 = p2 * p
1173 E2 = Fsum(x * y, x * z, y * z, -p2 * _3_0)
1174 E2p = E2 * p
1175 xyz = x * y * z
1176 return _Horner(S(_6_0), sqrt(a) * a * m, E2,
1177 Fsum(p3 * _4_0, xyz, E2p * _2_0),
1178 Fsum(p3 * _3_0, E2p, xyz * _2_0).fmul(p),
1179 xyz * p2, *over) # Fsum
1181# **) MIT License
1182#
1183# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1184#
1185# Permission is hereby granted, free of charge, to any person obtaining a
1186# copy of this software and associated documentation files (the "Software"),
1187# to deal in the Software without restriction, including without limitation
1188# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1189# and/or sell copies of the Software, and to permit persons to whom the
1190# Software is furnished to do so, subject to the following conditions:
1191#
1192# The above copyright notice and this permission notice shall be included
1193# in all copies or substantial portions of the Software.
1194#
1195# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1196# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1197# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1198# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1199# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1200# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1201# OTHER DEALINGS IN THE SOFTWARE.