1

2# -*- coding: utf-8 -*-

3

4 u'''I{Karney}'s elliptic functions and integrals.

5

6 Class L{Elliptic} transcoded from I{Charles Karney}'s C++ class U{EllipticFunction

7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1EllipticFunction.html>}

8 to pure Python, including symmetric integrals L{Elliptic.fRC}, L{Elliptic.fRD},

9 L{Elliptic.fRF}, L{Elliptic.fRG} and L{Elliptic.fRJ} as C{static methods}.

10

11 Python method names follow the C++ member functions, I{except}:

12

13 - member functions I{without arguments} are mapped to Python properties

14 prefixed with C{"c"}, for example C{E()} is property C{cE},

15

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)},

19

20 - other Python method names conventionally start with a lower-case

21 letter or an underscore if private.

22

23 Following is a copy of I{Karney}'s U{EllipticFunction.hpp

24<https://GeographicLib.SourceForge.io/C++/doc/EllipticFunction_8hpp_source.html>}

25 file C{Header}.

26

27 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2023)

28 and licensed under the MIT/X11 License. For more information, see the

29 U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.

30

31 B{Elliptic integrals and functions.}

32

33 This provides the elliptic functions and integrals needed for

34 C{Ellipsoid}, C{GeodesicExact}, and C{TransverseMercatorExact}. Two

35 categories of function are provided:

36

37 - functions to compute U{symmetric elliptic integrals

38 <https://DLMF.NIST.gov/19.16.i>}

39

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>}.

43

44 In the latter case, an object is constructed giving the modulus

45 C{k} (and optionally the parameter C{alpha}). The modulus (and

46 parameter) are always passed as squares which allows C{k} to be

47 pure imaginary. (Confusingly, Abramowitz and Stegun call C{m = k**2}

48 the "parameter" and C{n = alpha**2} the "characteristic".)

49

50 In geodesic applications, it is convenient to separate the incomplete

51 integrals into secular and periodic components, e.g.

52

53 I{C{E(phi, k) = (2 E(k) / pi) [ phi + delta E(phi, k) ]}}

54

55 where I{C{delta E(phi, k)}} is an odd periodic function with

56 period I{C{pi}}.

57

58 The computation of the elliptic integrals uses the algorithms given

59 in 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,

62 13--26 (1995) with the additional optimizations given U{here

63<https://DLMF.NIST.gov/19.36.i>}.

64

65 The computation of the Jacobi elliptic functions uses the algorithm

66 given in U{R. Bulirsch, Numerical Calculation of Elliptic Integrals

67 and Elliptic Functions<https://DOI.org/10.1007/BF01397975>},

68 Numerische Mathematik 7, 78--90 (1965).

69

70 The notation follows U{NIST Digital Library of Mathematical Functions

71<https://DLMF.NIST.gov>} chapters U{19<https://DLMF.NIST.gov/19>} and

72 U{22<https://DLMF.NIST.gov/22>}.

73'''

74# make sure int/int division yields float quotient, see .basics

75 from __future__ import division as _; del _ # PYCHOK semicolon

76

77 from pygeodesy.basics import copysign0, map2, neg

78 from 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

81 from pygeodesy.errors import _ValueError, _xkwds_pop

82 from pygeodesy.fmath import fdot, Fsum, hypot1

83# from pygeodesy.fsums import Fsum # from .fmath

84 from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _SPACE_

85 from pygeodesy.karney import _ALL_LAZY, _signBit

86# from pygeodesy.lazily import _ALL_LAZY # from .karney

87 from pygeodesy.named import _Named, _NamedTuple, _NotImplemented

88 from pygeodesy.props import _allPropertiesOf_n, Property_RO, property_RO, \

89 _update_all

90 from pygeodesy.streprs import Fmt, unstr

91 from pygeodesy.units import Scalar, Scalar_

92 from pygeodesy.utily import sincos2, sincos2d

93

94 from math import asinh, atan, atan2, ceil, cosh, fabs, floor, sin, sqrt, tanh

95

96__all__ = _ALL_LAZY.elliptic

97__version__ = '23.06.08'

98

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

104

105

106 class _Complete(object):

107 '''(INTERAL) Hold complete integrals.

108 '''

109 def __init__(self, **kwds):

110 self.__dict__ = kwds

111

112

113 class _Deferred(list):

114 '''(INTERNAL) Collector for L{_Deferred_Fsum}.

115 '''

116 def __init__(self, *xs):

117 if xs:

118 list.__init__(self, xs)

119

120 def __add__(self, other): # PYCHOK no cover

121 return _NotImplemented(self, other)

122

123 def __iadd__(self, x): # overide C{list} += C{list}

124 # assert isscalar(x)

125 self.append(float(x))

126 return self

127

128 def __imul__(self, other): # PYCHOK no cover

129 return _NotImplemented(self, other)

130

131 def __isub__(self, x): # PYCHOK no cover

132 # assert isscalar(x)

133 self.append(-float(x))

134 return self

135

136 def __mul__(self, other): # PYCHOK no cover

137 return _NotImplemented(self, other)

138

139 def __radd__(self, other): # PYCHOK no cover

140 return _NotImplemented(self, other)

141

142 def __rsub__(self, other): # PYCHOK no cover

143 return _NotImplemented(self, other)

144

145 def __sub__(self, other): # PYCHOK no cover

146 return _NotImplemented(self, other)

147

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

152

153

154 class _Deferred_Fsum(Fsum):

155 '''(INTERNAL) Deferred L{Fsum}.

156 '''

157 name = NN # pre-named, overridden below

158

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)

163

164_Deferred_Fsum.name = _Deferred_Fsum.__name__ # PYCHOK once

165

166

167 class Elliptic(_Named):

168 '''Elliptic integrals and functions.

169

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

178

179 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None, name=NN):

180 '''Constructor, specifying the C{modulus} and C{parameter}.

181

182 @kwarg name: Optional name (C{str}).

183

184 @see: Method L{Elliptic.reset} for further details.

185

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)

192

193 if name:

194 self.name = name

195

196 @Property_RO

197 def alpha2(self):

198 '''Get α^2, the square of the parameter (C{float}).

199 '''

200 return self._alpha2

201

202 @Property_RO

203 def alphap2(self):

204 '''Get α'^2, the square of the complementary parameter (C{float}).

205 '''

206 return self._alphap2

207

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

214

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

221

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

228

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

235

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

242

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

249

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

256

257 def deltaD(self, sn, cn, dn):

258 '''The periodic Jahnke's incomplete elliptic integral.

259

260 @arg sn: sin(φ).

261 @arg cn: cos(φ).

262 @arg dn: sqrt(1 − k2 * sin(2φ)).

263

264 @return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}).

265

266 @raise EllipticError: Invalid invokation or no convergence.

267 '''

268 return self._deltaX(sn, cn, dn, self.cD, self.fD)

269

270 def deltaE(self, sn, cn, dn):

271 '''The periodic incomplete integral of the second kind.

272

273 @arg sn: sin(φ).

274 @arg cn: cos(φ).

275 @arg dn: sqrt(1 − k2 * sin(2φ)).

276

277 @return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}).

278

279 @raise EllipticError: Invalid invokation or no convergence.

280 '''

281 return self._deltaX(sn, cn, dn, self.cE, self.fE)

282

283 def deltaEinv(self, stau, ctau):

284 '''The periodic inverse of the incomplete integral of the second kind.

285

286 @arg stau: sin(τ)

287 @arg ctau: cos(τ)

288

289 @return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}).

290

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

296

297 def deltaF(self, sn, cn, dn):

298 '''The periodic incomplete integral of the first kind.

299

300 @arg sn: sin(φ).

301 @arg cn: cos(φ).

302 @arg dn: sqrt(1 − k2 * sin(2φ)).

303

304 @return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}).

305

306 @raise EllipticError: Invalid invokation or no convergence.

307 '''

308 return self._deltaX(sn, cn, dn, self.cK, self.fF)

309

310 def deltaG(self, sn, cn, dn):

311 '''Legendre's periodic geodesic longitude integral.

312

313 @arg sn: sin(φ).

314 @arg cn: cos(φ).

315 @arg dn: sqrt(1 − k2 * sin(2φ)).

316

317 @return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}).

318

319 @raise EllipticError: Invalid invokation or no convergence.

320 '''

321 return self._deltaX(sn, cn, dn, self.cG, self.fG)

322

323 def deltaH(self, sn, cn, dn):

324 '''Cayley's periodic geodesic longitude difference integral.

325

326 @arg sn: sin(φ).

327 @arg cn: cos(φ).

328 @arg dn: sqrt(1 − k2 * sin(2φ)).

329

330 @return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}).

331

332 @raise EllipticError: Invalid invokation or no convergence.

333 '''

334 return self._deltaX(sn, cn, dn, self.cH, self.fH)

335

336 def deltaPi(self, sn, cn, dn):

337 '''The periodic incomplete integral of the third kind.

338

339 @arg sn: sin(φ).

340 @arg cn: cos(φ).

341 @arg dn: sqrt(1 − k2 * sin(2φ)).

342

343 @return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ

344 (C{float}).

345

346 @raise EllipticError: Invalid invokation or no convergence.

347 '''

348 return self._deltaX(sn, cn, dn, self.cPi, self.fPi)

349

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)

356

357 if _signBit(cn):

358 cn, sn = -cn, -sn

359 return fX(sn, cn, dn) * PI_2 / cX - atan2(sn, cn)

360

361 @Property_RO

362 def eps(self):

363 '''Get epsilon (C{float}).

364 '''

365 return self._reset_cDcEcKcKE_eps.eps

366

367 def fD(self, phi_or_sn, cn=None, dn=None):

368 '''Jahnke's incomplete elliptic integral in terms of

369 Jacobi elliptic functions.

370

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φ)).

374

375 @return: D(φ, k) as though φ ∈ (−π, π] (C{float}),

376 U{defined<https://DLMF.NIST.gov/19.2.E4>}.

377

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

385

386 return self._fXf(phi_or_sn, cn, dn, self.cD,

387 self.deltaD, _fD)

388

389 def fDelta(self, sn, cn):

390 '''The C{Delta} amplitude function.

391

392 @arg sn: sin(φ).

393 @arg cn: cos(φ).

394

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

401

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.

405

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φ)).

409

410 @return: E(φ, k) as though φ ∈ (−π, π] (C{float}),

411 U{defined<https://DLMF.NIST.gov/19.2.E5>}.

412

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

437

438 return self._fXf(phi_or_sn, cn, dn, self.cE,

439 self.deltaE, _fE)

440

441 def fEd(self, deg):

442 '''The incomplete integral of the second kind with

443 the argument given in degrees.

444

445 @arg deg: Angle (C{degrees}).

446

447 @return: E(π B{C{deg}}/180, k) (C{float}).

448

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

459

460 def fEinv(self, x):

461 '''The inverse of the incomplete integral of the second kind.

462

463 @arg x: Argument (C{float}).

464

465 @return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}}

466 (C{float}).

467

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

494

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.

498

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φ)).

502

503 @return: F(φ, k) as though φ ∈ (−π, π] (C{float}),

504 U{defined<https://DLMF.NIST.gov/19.2.E4>}.

505

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

513

514 return self._fXf(phi_or_sn, cn, dn, self.cK,

515 self.deltaF, _fF)

516

517 def fG(self, phi_or_sn, cn=None, dn=None):

518 '''Legendre's geodesic longitude integral in terms of

519 Jacobi elliptic functions.

520

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φ)).

524

525 @return: G(φ, k) as though φ ∈ (−π, π] (C{float}).

526

527 @raise EllipticError: Invalid invokation or no convergence.

528

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>}.

534

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)

541

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.

545

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φ)).

549

550 @return: H(φ, k) as though φ ∈ (−π, π] (C{float}).

551

552 @raise EllipticError: Invalid invokation or no convergence.

553

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>}.

558

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)

565

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.

569

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φ)).

573

574 @return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}).

575

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)

582

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

597

598 return self._fXf(phi_or_sn, cn, dn, cX, deltaX, _fX)

599

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)

613

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)

621

622 @Property_RO

623 def k2(self):

624 '''Get k^2, the square of the modulus (C{float}).

625 '''

626 return self._k2

627

628 @Property_RO

629 def kp2(self):

630 '''Get k'^2, the square of the complementary modulus (C{float}).

631 '''

632 return self._kp2

633

634 def reset(self, k2=0, alpha2=0, kp2=None, alphap2=None): # MCCABE 13

635 '''Reset the modulus, parameter and the complementaries.

636

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).

641

642 @raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}}

643 or B{C{alphap2}}.

644

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)

651

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)

654

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)

658

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)

675

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

705

706 return _Complete(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps)

707

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

744

745 return _Complete(cG=cG, cH=cH, cPi=cPi)

746

747 def sncndn(self, x):

748 '''The Jacobi elliptic function.

749

750 @arg x: The argument (C{float}).

751

752 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with

753 C{*n(B{x}, k)}.

754

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)

779

780 return Elliptic3Tuple(sn, cn, dn, iteration=self._iteration)

781

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)

792

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!

809

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))

815

816 @staticmethod

817 def fRC(x, y):

818 '''Degenerate symmetric integral of the first kind C{RC(x, y)}.

819

820 @return: C{RC(x, y)}, equivalent to C{RF(x, y, y)}.

821

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)

826

827 @staticmethod

828 def fRD(x, y, z):

829 '''Degenerate symmetric integral of the third kind C{RD(x, y, z)}.

830

831 @return: C{RD(x, y, z)}, equivalent to C{RJ(x, y, z, z)}.

832

833 @see: U{C{RD} definition<https://DLMF.NIST.gov/19.16.E5>} and

834 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.

835 '''

836 return _RD(None, x, y, z)

837

838 @staticmethod

839 def fRF(x, y, *z):

840 '''Symmetric or complete symmetric integral of the first kind

841 C{RF(x, y, z)} respectively C{RF(x, y)}.

842

843 @return: C{RF(x, y, z)} or C{RF(x, y)} for missing or zero B{C{z}}.

844

845 @see: U{C{RF} definition<https://DLMF.NIST.gov/19.16.E1>} and

846 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.

847 '''

848 return _RF3(None, x, y, *z) if z and z[0] else _RF2(None, x, y)

849

850 @staticmethod

851 def fRG(x, y, *z):

852 '''Symmetric or complete symmetric integral of the second kind

853 C{RG(x, y, z)} respectively C{RG(x, y)}.

854

855 @return: C{RG(x, y, z)} or C{RG(x, y)} for missing or zero B{C{z}}.

856

857 @see: U{C{RG} definition<https://DLMF.NIST.gov/19.16.E3>} and

858 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}.

859 '''

860 return _RG3(None, x, y, *z) if z and z[0] else (

861 _RG2(None, x, y) * _0_5)

862

863 @staticmethod

864 def fRJ(x, y, z, p):

865 '''Symmetric integral of the third kind C{RJ(x, y, z, p)}.

866

867 @return: C{RJ(x, y, z, p)}.

868

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)

873

874_allPropertiesOf_n(15, Elliptic) # # PYCHOK assert, see Elliptic.reset

875

876

877 class EllipticError(_ValueError):

878 '''Elliptic integral, function, convergence or other L{Elliptic} issue.

879 '''

880 pass

881

882

883 class Elliptic3Tuple(_NamedTuple):

884 '''3-Tuple C{(sn, cn, dn)} all C{scalar}.

885 '''

886 _Names_ = ('sn', 'cn', 'dn')

887 _Units_ = ( Scalar, Scalar, Scalar)

888

889

890 class _Lxyz(list):

891 '''(INTERNAL) Helper for C{_RD}, C{_RF3} and C{_RJ}.

892 '''

893 _a = None

894 _a0 = None

895

896 def __init__(self, *xyz_): # x, y, z [, p]

897 list.__init__(self, xyz_)

898

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 self._a0 = self._a = Fsum(*t).fover(n)

907 return self._a0

908

909 def asr3(self, a):

910 '''Compute next C{a}, C{sqrt(xyz_)} and C{fdot(sqrt(xyz))}.

911 '''

912 L = self

913 # assert a is L._a

914 s = map2(sqrt, L) # sqrt(x), srqt(y), sqrt(z) [, sqrt(p)]

915 r = fdot(s[:3], s[1], s[2], s[0]) # sqrt(x) * sqrt(y) + ...

916 L[:] = [(x + r) * _0_25 for x in L]

917 # assert L is self

918 L._a = a = (a + r) * _0_25

919 return a, s, r

920

921 def rescale(self, am, *xy_):

922 '''Rescale C{x}, C{y}, ...

923 '''

924 for x in xy_:

925 yield (self._a0 - x) / am

926

927 def thresh(self, Tol):

928 '''Return the convergence threshold.

929 '''

930 return max(fabs(self._a0 - x) for x in self) / Tol

931

932

933 def _horner(S, e1, E2, E3, E4, E5, *over):

934 '''(INTERNAL) Horner form for C{_RD} and C{_RJ} below.

935 '''

936 E22 = E2**2

937 # Polynomial is <https://DLMF.NIST.gov/19.36.E2>

938 # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52

939 # + 3*E5/26 - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20

940 # + 45*E2**2*E3/272 - 9*(E3*E4+E2*E5)/68)

941 # converted to Horner-like form ...

942 e = e1 * 4084080

943 S *= e

944 S += Fsum( E2 * -540540, 471240).fmul(E5)

945 S += Fsum( E3 * -540540, E2 * 612612, -556920).fmul(E4)

946 S += Fsum(E22 * 675675, E3 * 306306, E2 * -706860, 680680).fmul(E3)

947 S += Fsum(E22 * -255255, E2 * 417690, -875160).fmul(E2)

948 S += 4084080

949 return S.fover((e * over[0]) if over else e)

950

951

952 def _invokationError(name, *args): # PYCHOK no cover

953 '''(INTERNAL) Return an L{EllipticError}.

954 '''

955 n = _DOT_(Elliptic.__name__, name)

956 n = _SPACE_(_invokation_, n)

957 return EllipticError(NN(n, repr(args))) # unstr

958

959

960 def _iterations(inst, i):

961 '''(INTERNAL) Aggregate iterations B{C{i}}.

962 '''

963 if inst:

964 inst._iteration += i

965

966

967 def _no_convergenceError(d, tol, where, *args, **kwds_thresh): # PYCHOK no cover

968 '''(INTERNAL) Return an L{EllipticError}.

969 '''

970 n = Elliptic.__name__

971 if where:

972 n = _DOT_(n, where.__name__)

973 if kwds_thresh:

974 q = _xkwds_pop(kwds_thresh, thresh=False)

975 t = unstr(n, *args, **kwds_thresh)

976 else:

977 q = False

978 t = unstr(n, *args)

979 return EllipticError(Fmt.no_convergence(d, tol, thresh=q), txt=t)

980

981

982 def _3over(a, b):

983 '''(INTERNAL) Return C{3 / (a * b)}.

984 '''

985 return _3_0 / (a * b)

986

987

988 def _RC(unused, x, y):

989 '''(INTERNAL) Defined only for y != 0 and x >= 0.

990 '''

991 d = x - y

992 if d < 0: # catch _NaN

993 # <https://DLMF.NIST.gov/19.2.E18>

994 d = -d

995 r = atan(sqrt(d / x)) if x > 0 else PI_2

996 elif d == _0_0: # XXX d < EPS0? or EPS02 or _EPSmin

997 d, r = y, _1_0

998 elif y > 0: # <https://DLMF.NIST.gov/19.2.E19>

999 r = asinh(sqrt(d / y)) # atanh(sqrt((x - y) / x))

1000 elif y < 0: # <https://DLMF.NIST.gov/19.2.E20>

1001 r = asinh(sqrt(-x / y)) # atanh(sqrt(x / (x - y)))

1002 else:

1003 raise _invokationError(Elliptic.fRC.__name__, x, y)

1004 return r / sqrt(d)

1005

1006

1007 def _RD(inst, x, y, z, *over):

1008 '''(INTERNAL) Carlson, eqs 2.28 - 2.34.

1009 '''

1010 L = _Lxyz(x, y, z)

1011 a = L.a0(5)

1012 q = L.thresh(_TolRF)

1013 S = _Deferred()

1014 am, m = a, 1

1015 for i in range(_TRIPS):

1016 if fabs(am) > q: # max 7 trips

1017 _iterations(inst, i)

1018 break

1019 t = L[2] # z0...n

1020 a, s, r = L.asr3(a)

1021 S += _1_0 / ((t + r) * s[2] * m)

1022 m *= 4

1023 am = a * m

1024 else: # PYCHOK no cover

1025 raise _no_convergenceError(am, q, Elliptic.fRD, x, y, z, *over,

1026 thresh=True)

1027 x, y = L.rescale(-am, x, y)

1028 xy = x * y

1029 z = (x + y) / _3_0

1030 z2 = z**2

1031 S = S.Fsum.fmul(_3_0)

1032 return _horner(S, am * sqrt(a),

1033 xy - _6_0 * z2,

1034 (xy * _3_0 - _8_0 * z2) * z,

1035 (xy - z2) * _3_0 * z2,

1036 xy * z2 * z, *over)

1037

1038

1039 def _RF2(inst, x, y): # 2-arg version, z=0

1040 '''(INTERNAL) Carlson, eqs 2.36 - 2.38.

1041 '''

1042 a, b = sqrt(x), sqrt(y)

1043 if a < b:

1044 a, b = b, a

1045 for i in range(_TRIPS):

1046 t = _TolRG0 * a

1047 if fabs(a - b) <= t: # max 4 trips

1048 _iterations(inst, i)

1049 return (PI / (a + b))

1050 a, b = ((a + b) * _0_5), sqrt(a * b)

1051 else: # PYCHOK no cover

1052 raise _no_convergenceError(a - b, t, Elliptic.fRF, x, y)

1053

1054

1055 def _RF3(inst, x, y, z): # 3-arg version

1056 '''(INTERNAL) Carlson, eqs 2.2 - 2.7.

1057 '''

1058 L = _Lxyz(x, y, z)

1059 a = L.a0(3)

1060 q = L.thresh(_TolRF)

1061 am, m = a, 1

1062 for i in range(_TRIPS):

1063 if fabs(am) > q: # max 6 trips

1064 _iterations(inst, i)

1065 break

1066 a, _, _ = L.asr3(a)

1067 m *= 4

1068 am = a * m

1069 else: # PYCHOK no cover

1070 raise _no_convergenceError(am, q, Elliptic.fRF, x, y, z,

1071 thresh=True)

1072 x, y = L.rescale(am, x, y)

1073 z = neg(x + y)

1074 xy = x * y

1075 e2 = xy - z**2

1076 e3 = xy * z

1077 e4 = e2**2

1078 # Polynomial is <https://DLMF.NIST.gov/19.36.E1>

1079 # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44

1080 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16)

1081 # converted to Horner-like form ...

1082 S = Fsum(e4 * 15015, e3 * 6930, e2 * -16380, 17160).fmul(e3)

1083 S += Fsum(e4 * -5775, e2 * 10010, -24024).fmul(e2)

1084 S += Fsum(240240)

1085 return S.fover(sqrt(a) * 240240)

1086

1087

1088 def _RG2(inst, x, y): # 2-args and I{doubled}

1089 '''(INTERNAL) Carlson, eqs 2.36 - 2.39.

1090 '''

1091 a, b = sqrt(x), sqrt(y)

1092 if a < b:

1093 a, b = b, a

1094 ab = a - b # fabs(a - b)

1095 S = _Deferred(_0_5 * (a + b)**2)

1096 m = -1

1097 for i in range(_TRIPS): # max 4 trips

1098 t = _TolRG0 * a

1099 if ab <= t:

1100 _iterations(inst, i)

1101 return S.Fsum.fover((a + b) / PI_2)

1102 a, b = ((a + b) * _0_5), sqrt(a * b)

1103 ab = fabs(a - b)

1104 S += ab**2 * m

1105 m *= 2

1106 else: # PYCHOK no cover

1107 raise _no_convergenceError(ab, t, Elliptic.fRG, x, y)

1108

1109

1110 def _RG3(inst, x, y, z): # 3-arg version

1111 '''(INTERNAL) Never called with zero B{C{z}}, see C{.fRG}.

1112 '''

1113# if not z:

1114# y, z = z, y

1115 rd = (x - z) * (z - y) # - (y - z)

1116 if rd: # Carlson, eq 1.7

1117 rd = _RD(inst, x, y, z, _3_0 * z / rd)

1118 xyz = x * y

1119 if xyz:

1120 xyz = sqrt(xyz / z**3)

1121 return Fsum(_RF3(inst, x, y, z), rd, xyz).fover(_2_0 / z)

1122

1123

1124 def _RJ(inst, x, y, z, p, *over):

1125 '''(INTERNAL) Carlson, eqs 2.17 - 2.25.

1126 '''

1127 def _xyzp(x, y, z, p):

1128 return (x + p) * (y + p) * (z + p)

1129

1130 L = _Lxyz(x, y, z, p)

1131 a = L.a0(5)

1132 q = L.thresh(_TolRD)

1133 S = _Deferred()

1134 n = neg(_xyzp(x, y, z, -p))

1135 am, m = a, 1

1136 for i in range(_TRIPS):

1137 if fabs(am) > q: # max 7 trips

1138 _iterations(inst, i)

1139 break

1140 a, s, _ = L.asr3(a)

1141 d = _xyzp(*s)

1142 if n:

1143 rc = _RC(inst, _1_0, n / d**2 + _1_0)

1144 n *= _1_64th # /= chokes PyChecker

1145 else:

1146 rc = _1_0 # == _RC(None, _1_0, _1_0)

1147 S += rc / (d * m)

1148 m *= 4

1149 am = a * m

1150 else: # PYCHOK no cover

1151 raise _no_convergenceError(am, q, Elliptic.fRJ, x, y, z, p,

1152 thresh=True)

1153 x, y, z = L.rescale(am, x, y, z)

1154 xyz = x * y * z

1155 p = Fsum(x, y, z).fover(_N_2_0)

1156 p2 = p**2

1157 p3 = p2*p

1158 E2 = Fsum(x * y, x * z, y * z, -p2 * _3_0)

1159 E2p = E2 * p

1160 S = S.Fsum.fmul(_6_0)

1161 return _horner(S, am * sqrt(a), E2,

1162 Fsum(p3 * _4_0, xyz, E2p * _2_0),

1163 Fsum(p3 * _3_0, E2p, xyz * _2_0).fmul(p),

1164 xyz * p2, *over)

1165

1166# **) MIT License

1167#

1168# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

1169#

1170# Permission is hereby granted, free of charge, to any person obtaining a

1171# copy of this software and associated documentation files (the "Software"),

1172# to deal in the Software without restriction, including without limitation

1173# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1174# and/or sell copies of the Software, and to permit persons to whom the

1175# Software is furnished to do so, subject to the following conditions:

1176#

1177# The above copyright notice and this permission notice shall be included

1178# in all copies or substantial portions of the Software.

1179#

1180# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1181# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1182# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1183# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1184# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1185# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1186# OTHER DEALINGS IN THE SOFTWARE.