1 2# -*- coding: utf-8 -*- 3 4u'''I{Karney}'s elliptic functions and integrals. 5 6Class L{Elliptic} transcoded from I{Charles Karney}'s C++ class U{EllipticFunction 7<https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1EllipticFunction.html>} 8to pure Python. 9 10Python method names follow the C++ member functions, I{except}: 11 12 - member functions I{without arguments} are mapped to Python properties 13 prefixed with C{"c"}, for example C{E()} is property C{cE}, 14 15 - member functions with 1 or 3 arguments are renamed to Python methods 16 starting with an C{"f"}, example C{E(psi)} to C{fE(psi)} and C{E(sn, 17 cn, dn)} to C{fE(sn, cn, dn)}, 18 19 - other Python method names conventionally start with a lower-case 20 letter or an underscore if private. 21 22Following is a copy of I{Karney}'s U{EllipticFunction.hpp 23<https://GeographicLib.SourceForge.io/html/EllipticFunction_8hpp_source.html>} 24file C{Header}. 25 26Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2021) 27and licensed under the MIT/X11 License. For more information, see the 28U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 29 30B{Elliptic integrals and functions.} 31 32This provides the elliptic functions and integrals needed for 33C{Ellipsoid}, C{GeodesicExact}, and C{TransverseMercatorExact}. Two 34categories of function are provided: 35 36 - functions to compute U{symmetric elliptic integrals 37 <https://DLMF.NIST.gov/19.16.i>} 38 39 - methods to compute U{Legrendre's elliptic integrals 40 <https://DLMF.NIST.gov/19.2.ii>} and U{Jacobi elliptic 41 functions<https://DLMF.NIST.gov/22.2>}. 42 43In the latter case, an object is constructed giving the modulus 44C{k} (and optionally the parameter C{alpha}). The modulus (and 45parameter) are always passed as squares which allows C{k} to be 46pure imaginary. (Confusingly, Abramowitz and Stegun call C{m = k**2} 47the "parameter" and C{n = alpha**2} the "characteristic".) 48 49In geodesic applications, it is convenient to separate the incomplete 50integrals into secular and periodic components, e.g. 51 52I{C{E(phi, k) = (2 E(k) / pi) [ phi + delta E(phi, k) ]}} 53 54where I{C{delta E(phi, k)}} is an odd periodic function with 55period I{C{pi}}. 56 57The computation of the elliptic integrals uses the algorithms given 58in U{B. C. Carlson, Computation of real or complex elliptic integrals 59<https://DOI.org/10.1007/BF02198293>} (also available U{here 60<https://ArXiv.org/pdf/math/9409227.pdf>}), Numerical Algorithms 10, 6113--26 (1995) with the additional optimizations given U{here 62<https://DLMF.NIST.gov/19.36.i>}. 63 64The computation of the Jacobi elliptic functions uses the algorithm 65given in U{R. Bulirsch, Numerical Calculation of Elliptic Integrals 66and Elliptic Functions<https://DOI.org/10.1007/BF01397975>}, 67Numerische Mathematik 7, 78--90 (1965). 68 69The notation follows U{NIST Digital Library of Mathematical Functions 70<https://DLMF.NIST.gov>} chapters U{19<https://DLMF.NIST.gov/19>} and 71U{22<https://DLMF.NIST.gov/22>}. 72''' 73# make sure int/int division yields float quotient, see .basics 74from __future__ import division 75 76from pygeodesy.basics import copysign0, map2, neg 77from pygeodesy.errors import _ValueError 78from pygeodesy.fmath import fdot, fmean_, Fsum, fsum_, hypot1 79from pygeodesy.interns import EPS, INF, NN, PI, PI_2, PI_4, \ 80 _EPStol as _TolJAC, _convergence_, \ 81 _no_, _SPACE_, _0_0, _0_125, _0_25, \ 82 _0_5, _1_0, _2_0, _3_0, _4_0, _5_0, \ 83 _6_0, _8_0, _180_0, _360_0 84from pygeodesy.lazily import _ALL_LAZY 85from pygeodesy.named import _Named, _NamedTuple 86from pygeodesy.props import Property_RO, property_RO, _update_all 87# from pygeodesy.streprs import unstr 88from pygeodesy.units import Scalar, Scalar_ 89from pygeodesy.utily import sincos2, sincos2d 90 91from math import asinh, atan, atan2, ceil, cosh, floor, sin, \ 92 sqrt, tanh 93 94__all__ = _ALL_LAZY.elliptic 95__version__ = '21.07.31' 96 97_TolRD = pow(EPS * 0.002, _0_125) 98_TolRF = pow(EPS * 0.030, _0_125) 99_TolRG0 = _TolJAC * 2.7 100_TRIPS = 15 # Max depth, 7 might be enough, by .etm 101 102 103class EllipticError(_ValueError): 104 '''Elliptic integral, function, convergence or other L{Elliptic} issue. 105 ''' 106 pass 107 108 109class Elliptic3Tuple(_NamedTuple): 110 '''3-Tuple C{(sn, cn, dn)} all C{scalar}. 111 ''' 112 _Names_ = ('sn', 'cn', 'dn') 113 _Units_ = ( Scalar, Scalar, Scalar) 114 115 116class Elliptic(_Named): 117 '''Elliptic integrals and functions. 118 119 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/ 120 html/classGeographicLib_1_1EllipticFunction.html#details>}. 121 ''' 122 _alpha2 = 0 123 _alphap2 = 0 124 _eps = EPS 125 _iteration = None # only .fEinv and .sncndn 126 _k2 = 0 127 _kp2 = 0 128 129 _cD = INF 130 _cE = _1_0 131 _cG = INF 132 _cH = INF 133 _cK = INF 134 _cKE = INF 135 _cPi = INF 136 137 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None): 138 '''Constructor, specifying the C{modulus} and C{parameter}. 139 140 @kwarg k2: Modulus squared (C{float}, 0 <= k^2 <= 1). 141 @kwarg alpha2: Parameter squared (C{float}, 0 <= α^2 <= 1). 142 @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0). 143 @kwarg alphap2: Complementary parameter squared (C{float}, α'^2 >= 0). 144 145 @see: Method L{reset} for further details. 146 147 @note: If only elliptic integrals of the first and second kinds 148 are needed, use C{B{alpha2}=0}, the default value. In 149 that case, we have C{Π(φ, 0, k) = F(φ, k), G(φ, 0, k) = 150 E(φ, k)} and C{H(φ, 0, k) = F(φ, k) - D(φ, k)}. 151 ''' 152 self.reset(k2=k2, alpha2=alpha2, kp2=kp2, alphap2=alphap2) 153 154 @Property_RO 155 def alpha2(self): 156 '''Get α^2, the square of the parameter (C{float}). 157 ''' 158 return self._alpha2 159 160 @Property_RO 161 def alphap2(self): 162 '''Get α'^2, the square of the complementary parameter (C{float}). 163 ''' 164 return self._alphap2 165 166 @Property_RO 167 def cD(self): 168 '''Get Jahnke's complete integral C{D(k)} (C{float}), 169 U{defined<https://DLMF.NIST.gov/19.2.E6>}. 170 ''' 171 return self._cD 172 173 @Property_RO 174 def cE(self): 175 '''Get the complete integral of the second kind C{E(k)} 176 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}. 177 ''' 178 return self._cE 179 180 @Property_RO 181 def cG(self): 182 '''Get Legendre's complete geodesic longitude integral 183 C{G(α^2, k)} (C{float}). 184 ''' 185 return self._cG 186 187 @Property_RO 188 def cH(self): 189 '''Get Cayley's complete geodesic longitude difference integral 190 C{H(α^2, k)} (C{float}). 191 ''' 192 return self._cH 193 194 @Property_RO 195 def cK(self): 196 '''Get the complete integral of the first kind C{K(k)} 197 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}. 198 ''' 199 return self._cK 200 201 @Property_RO 202 def cKE(self): 203 '''Get the difference between the complete integrals of the 204 first and second kinds, C{K(k) − E(k)} (C{float}). 205 ''' 206 return self._cKE 207 208 @Property_RO 209 def cPi(self): 210 '''Get the complete integral of the third kind C{Pi(α^2, k)} 211 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}. 212 ''' 213 return self._cPi 214 215 def deltaD(self, sn, cn, dn): 216 '''The periodic Jahnke's incomplete elliptic integral. 217 218 @arg sn: sin(φ). 219 @arg cn: cos(φ). 220 @arg dn: sqrt(1 − k2 sin(2φ)). 221 222 @return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}). 223 224 @raise EllipticError: Invalid invokation or no convergence. 225 ''' 226 return self._deltaX(sn, cn, dn, self.cD, self.fD) 227 228 def deltaE(self, sn, cn, dn): 229 '''The periodic incomplete integral of the second kind. 230 231 @arg sn: sin(φ). 232 @arg cn: cos(φ). 233 @arg dn: sqrt(1 − k2 sin(2φ)). 234 235 @return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}). 236 237 @raise EllipticError: Invalid invokation or no convergence. 238 ''' 239 return self._deltaX(sn, cn, dn, self.cE, self.fE) 240 241 def deltaEinv(self, stau, ctau): 242 '''The periodic inverse of the incomplete integral of the second kind. 243 244 @arg stau: sin(τ) 245 @arg ctau: cos(τ) 246 247 @return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}). 248 249 @raise EllipticError: No convergence. 250 ''' 251 # Function is periodic with period pi 252 t = atan2(-stau, -ctau) if ctau < 0 else atan2(stau, ctau) 253 return self.fEinv(t * self._cE / PI_2) - t 254 255 def deltaF(self, sn, cn, dn): 256 '''The periodic incomplete integral of the first kind. 257 258 @arg sn: sin(φ). 259 @arg cn: cos(φ). 260 @arg dn: sqrt(1 − k2 sin(2φ)). 261 262 @return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}). 263 264 @raise EllipticError: Invalid invokation or no convergence. 265 ''' 266 return self._deltaX(sn, cn, dn, self.cK, self.fF) 267 268 def deltaG(self, sn, cn, dn): 269 '''Legendre's periodic geodesic longitude integral. 270 271 @arg sn: sin(φ). 272 @arg cn: cos(φ). 273 @arg dn: sqrt(1 − k2 sin(2φ)). 274 275 @return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}). 276 277 @raise EllipticError: Invalid invokation or no convergence. 278 ''' 279 return self._deltaX(sn, cn, dn, self.cG, self.fG) 280 281 def deltaH(self, sn, cn, dn): 282 '''Cayley's periodic geodesic longitude difference integral. 283 284 @arg sn: sin(φ). 285 @arg cn: cos(φ). 286 @arg dn: sqrt(1 − k2 sin(2φ)). 287 288 @return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}). 289 290 @raise EllipticError: Invalid invokation or no convergence. 291 ''' 292 return self._deltaX(sn, cn, dn, self.cH, self.fH) 293 294 def deltaPi(self, sn, cn, dn): 295 '''The periodic incomplete integral of the third kind. 296 297 @arg sn: sin(φ). 298 @arg cn: cos(φ). 299 @arg dn: sqrt(1 − k2 sin(2φ)). 300 301 @return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ 302 (C{float}). 303 304 @raise EllipticError: Invalid invokation or no convergence. 305 ''' 306 return self._deltaX(sn, cn, dn, self.cPi, self.fPi) 307 308 def _deltaX(self, sn, cn, dn, cX, fX): 309 '''(INTERNAL) Helper for C{.deltaD} thru C{.deltaPi}. 310 ''' 311 if None in (cn, dn): 312 t = self.classname + '.delta' + fX.__name__[1:] 313 raise _invokationError(t, sn, cn, dn) 314 315 if cn < 0: 316 cn, sn = -cn, -sn 317 return fX(sn, cn, dn) * PI_2 / cX - atan2(sn, cn) 318 319 @Property_RO 320 def eps(self): 321 '''Get epsilon (C{float}). 322 ''' 323 return self._eps 324 325 def fD(self, phi_or_sn, cn=None, dn=None): 326 '''Jahnke's incomplete elliptic integral in terms of 327 Jacobi elliptic functions. 328 329 @arg phi_or_sn: φ or sin(φ). 330 @kwarg cn: C{None} or cos(φ). 331 @kwarg dn: C{None} or sqrt(1 − k2 sin(2φ)). 332 333 @return: D(φ, k) as though φ ∈ (−π, π] (C{float}), 334 U{defined<https://DLMF.NIST.gov/19.2.E4>}. 335 336 @raise EllipticError: Invalid invokation or no convergence. 337 ''' 338 def _fD(sn, cn, dn): 339 return abs(sn) * sn**2 * _RD_3(cn**2, dn**2, _1_0) 340 341 return self._fXf(phi_or_sn, cn, dn, self.cD, 342 self.deltaD, _fD) 343 344 def fDelta(self, sn, cn): 345 '''The C{Delta} amplitude function. 346 347 @arg sn: sin(φ). 348 @arg cn: cos(φ). 349 350 @return: sqrt(1 − k2 sin(2φ)) (C{float}). 351 ''' 352 k2 = self.k2 353 return sqrt((_1_0 - k2 * sn**2) if k2 < 0 else 354 (self.kp2 + k2 * cn**2)) 355 356 def fE(self, phi_or_sn, cn=None, dn=None): 357 '''The incomplete integral of the second kind in terms of 358 Jacobi elliptic functions. 359 360 @arg phi_or_sn: φ or sin(φ). 361 @kwarg cn: C{None} or cos(φ). 362 @kwarg dn: C{None} or sqrt(1 − k2 sin(2φ)). 363 364 @return: E(φ, k) as though φ ∈ (−π, π] (C{float}), 365 U{defined<https://DLMF.NIST.gov/19.2.E5>}. 366 367 @raise EllipticError: Invalid invokation or no convergence. 368 ''' 369 def _fE(sn, cn, dn): 370 sn2, cn2, dn2 = sn**2, cn**2, dn**2 371 kp2, k2 = self.kp2, self.k2 372 if k2 <= 0: # Carlson, eq. 4.6, <https://DLMF.NIST.gov/19.25.E9> 373 ei = _RF(cn2, dn2, _1_0) - k2 * sn2 * _RD_3(cn2, dn2, _1_0) 374 elif kp2 >= 0: # <https://DLMF.NIST.gov/19.25.E10> 375 ei = fsum_(kp2 * _RF(cn2, dn2, _1_0), 376 kp2 * k2 * sn2 * _RD_3(cn2, _1_0, dn2), 377 k2 * abs(cn) / dn) 378 else: # <https://DLMF.NIST.gov/19.25.E11> 379 ei = dn / abs(cn) - kp2 * sn2 * _RD_3(dn2, _1_0, cn2) 380 return ei * abs(sn) 381 382 return self._fXf(phi_or_sn, cn, dn, self.cE, 383 self.deltaE, _fE) 384 385 def fEd(self, deg): 386 '''The incomplete integral of the second kind with 387 the argument given in degrees. 388 389 @arg deg: Angle (C{degrees}). 390 391 @return: E(π B{C{deg}}/180, k) (C{float}). 392 393 @raise EllipticError: No convergence. 394 ''' 395 if abs(deg) < _180_0: 396 e = _0_0 397 else: 398 n = ceil(deg / _360_0 - _0_5) 399 deg -= n * _360_0 400 e = n * _4_0 * self.cE 401 sn, cn = sincos2d(deg) 402 return self.fE(sn, cn, self.fDelta(sn, cn)) + e 403 404 def fEinv(self, x): 405 '''The inverse of the incomplete integral of the second kind. 406 407 @arg x: Argument (C{float}). 408 409 @return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}} 410 (C{float}). 411 412 @raise EllipticError: No convergence. 413 ''' 414 E2 = self.cE * _2_0 415 n = floor(x / E2 + _0_5) 416 x -= E2 * n # x now in [-ec, ec) 417 # linear approximation 418 phi = PI * x / E2 # phi in [-pi/2, pi/2) 419 Phi = Fsum(phi) 420 # first order correction 421 phi = Phi.fsum_(-self.eps * sin(2 * phi) * _0_5) 422 # For kp2 close to zero use asin(x/.cE) or J. P. Boyd, 423 # Applied Math. and Computation 218, 7005-7013 (2012) 424 # <https://DOI.org/10.1016/j.amc.2011.12.021> 425 for self._iteration in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC 426 sn, cn, dn = self._sncndn3(phi) 427 phi, e = Phi.fsum2_((x - self.fE(sn, cn, dn)) / dn) 428 if abs(e) < _TolJAC: 429 if n: 430 phi = Phi.fsum_(n * PI) 431 return phi 432 raise _convergenceError(self.fEinv, x) 433 434 def fF(self, phi_or_sn, cn=None, dn=None): 435 '''The incomplete integral of the first kind in terms of 436 Jacobi elliptic functions. 437 438 @arg phi_or_sn: φ or sin(φ). 439 @kwarg cn: C{None} or cos(φ). 440 @kwarg dn: C{None} or sqrt(1 − k2 sin(2φ)). 441 442 @return: F(φ, k) as though φ ∈ (−π, π] (C{float}), 443 U{defined<https://DLMF.NIST.gov/19.2.E4>}. 444 445 @raise EllipticError: Invalid invokation or no convergence. 446 ''' 447 def _fF(sn, cn, dn): 448 return abs(sn) * _RF(cn**2, dn**2, _1_0) 449 450 return self._fXf(phi_or_sn, cn, dn, self.cK, 451 self.deltaF, _fF) 452 453 def fG(self, phi_or_sn, cn=None, dn=None): 454 '''Legendre's geodesic longitude integral in terms of 455 Jacobi elliptic functions. 456 457 @arg phi_or_sn: φ or sin(φ). 458 @kwarg cn: C{None} or cos(φ). 459 @kwarg dn: C{None} or sqrt(1 − k2 sin(2φ)). 460 461 @return: G(φ, k) as though φ ∈ (−π, π] (C{float}). 462 463 @raise EllipticError: Invalid invokation or no convergence. 464 465 @note: Legendre expresses the longitude of a point on 466 the geodesic in terms of this combination of 467 elliptic integrals in U{Exercices de Calcul 468 Intégral, Vol 1 (1811), p 181<https://Books. 469 Google.com/books?id=riIOAAAAQAAJ&pg=PA181>}. 470 471 @see: U{Geodesics in terms of elliptic integrals<https:// 472 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>} 473 for the expression for the longitude in terms of this function. 474 ''' 475 return self._fXa(phi_or_sn, cn, dn, self.alpha2 - self.k2, 476 self.cG, self.deltaG) 477 478 def fH(self, phi_or_sn, cn=None, dn=None): 479 '''Cayley's geodesic longitude difference integral in terms of 480 Jacobi elliptic functions. 481 482 @arg phi_or_sn: φ or sin(φ). 483 @kwarg cn: C{None} or cos(φ). 484 @kwarg dn: C{None} or sqrt(1 − k2 sin(2φ)). 485 486 @return: H(φ, k) as though φ ∈ (−π, π] (C{float}). 487 488 @raise EllipticError: Invalid invokation or no convergence. 489 490 @note: Cayley expresses the longitude difference of a point 491 on the geodesic in terms of this combination of 492 elliptic integrals in U{Phil. Mag. B{40} (1870), p 333 493 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}. 494 495 @see: U{Geodesics in terms of elliptic integrals<https:// 496 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>} 497 for the expression for the longitude in terms of this function. 498 ''' 499 return self._fXa(phi_or_sn, cn, dn, -self.alphap2, 500 self.cH, self.deltaH) 501 502 def fPi(self, phi_or_sn, cn=None, dn=None): 503 '''The incomplete integral of the third kind in terms of 504 Jacobi elliptic functions. 505 506 @arg phi_or_sn: φ or sin(φ). 507 @kwarg cn: C{None} or cos(φ). 508 @kwarg dn: C{None} or sqrt(1 − k2 sin(2φ)). 509 510 @return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}). 511 512 @raise EllipticError: Invalid invokation or no convergence. 513 ''' 514 return self._fXa(phi_or_sn, cn, dn, self.alpha2, 515 self.cPi, self.deltaPi) 516 517 def _fXa(self, phi_or_sn, cn, dn, aX, cX, deltaX): 518 '''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}. 519 ''' 520 def _fX(sn, cn, dn): 521 cn2, sn2, dn2 = cn**2, sn**2, dn**2 522 return abs(sn) * (_RF(cn2, dn2, _1_0) + aX * sn2 * 523 _RJ_3(cn2, dn2, _1_0, cn2 + self.alphap2 * sn2)) 524 525 return self._fXf(phi_or_sn, cn, dn, cX, deltaX, _fX) 526 527 def _fXf(self, phi_or_sn, cn, dn, cX, deltaX, fX): 528 '''(INTERNAL) Helper for C{f.D}, C{.fE}, C{.fF} and C{._fXa}. 529 ''' 530 phi = sn = phi_or_sn 531 if cn is dn is None: # fX(phi) call 532 sn, cn, dn = self._sncndn3(phi) 533 if abs(phi) >= PI: 534 return (deltaX(sn, cn, dn) + phi) * cX / PI_2 535 # fall through 536 elif None in (cn, dn): 537 t = self.classname + '.f' + deltaX.__name__[5:] 538 raise _invokationError(t, sn, cn, dn) 539 540 if cn < 0: # enforce usual trig-like symmetries 541 xi = _2_0 * cX - fX(sn, cn, dn) 542 elif cn > 0: 543 xi = fX(sn, cn, dn) 544 else: 545 xi = cX 546 return copysign0(xi, sn) 547 548 @property_RO 549 def iteration(self): 550 '''Get the most recent C{Elliptic.fEinv} or C{Elliptic.sncndn} 551 iteration number (C{int}) or C{None} if not available/applicable. 552 ''' 553 return self._iteration 554 555 @Property_RO 556 def k2(self): 557 '''Get k^2, the square of the modulus (C{float}). 558 ''' 559 return self._k2 560 561 @Property_RO 562 def kp2(self): 563 '''Get k'^2, the square of the complementary modulus (C{float}). 564 ''' 565 return self._kp2 566 567 def reset(self, k2=0, alpha2=0, kp2=None, alphap2=None): # MCCABE 13 568 '''Reset the modulus and parameter. 569 570 @kwarg k2: modulus squared (C{float}, -INF <= k^2<= 1). 571 @kwarg alpha2: parameter (C{float}, -INF <= α^2 <= 1). 572 @kwarg kp2: complementary modulus squared (C{float}, k'^2 >= 0). 573 @kwarg alphap2: complementary parameter squared (C{float}, α'^2 >= 0). 574 575 @raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}} 576 or B{C{alphap2}} or no convergence. 577 578 @note: The arguments must satisfy C{B{k2} + B{kp2} = 1} and 579 C{B{alpha2} + B{alphap2} = 1}. No checking is done 580 that these conditions are met to enable accuracy to 581 be maintained, e.g., when C{k} is very close to unity. 582 ''' 583 _update_all(self) 584 585 self._k2 = k2 = Scalar_(k2=k2, Error=EllipticError, low=None, high=_1_0) 586 587 self._alpha2 = alpha2 = Scalar_(alpha2=alpha2, Error=EllipticError, low=None, high=_1_0) 588 589 self._kp2 = kp2 = Scalar_(kp2=((_1_0 - k2) if kp2 is None else kp2), Error=EllipticError) 590 591 self._alphap2 = alphap2 = Scalar_(alphap2=((_1_0 - alpha2) if alphap2 is None else alphap2), 592 Error=EllipticError) 593 594 # Values of complete elliptic integrals for k = 0,1 and alpha = 0,1 595 # K E D 596 # k = 0: pi/2 pi/2 pi/4 597 # k = 1: inf 1 inf 598 # Pi G H 599 # k = 0, alpha = 0: pi/2 pi/2 pi/4 600 # k = 1, alpha = 0: inf 1 1 601 # k = 0, alpha = 1: inf inf pi/2 602 # k = 1, alpha = 1: inf inf inf 603 # 604 # G(0, k) = Pi(0, k) = H(1, k) = E(k) 605 # H(0, k) = K(k) - D(k) 606 # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2)) 607 # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1)) 608 # Pi(alpha2, 1) = inf 609 # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2) 610 if k2: 611 if kp2: 612 self._eps = k2 / (sqrt(kp2) + _1_0)**2 613 # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3 614 # <https://DLMF.NIST.gov/19.25.E1> 615 self._cD = _RD_3(_0_0, kp2, _1_0) 616 # Complete elliptic integral E(k), Carlson eq. 4.2 617 # <https://DLMF.NIST.gov/19.25.E1> 618 self._cE = _RG_(kp2, _1_0) * _2_0 619 # Complete elliptic integral K(k), Carlson eq. 4.1 620 # <https://DLMF.NIST.gov/19.25.E1> 621 self._cK = _RF_(kp2, _1_0) 622 self._cKE = k2 * self.cD 623 else: 624 self._eps = k2 625 self._cD = self._cK = self._cKE = INF 626 self._cE = _1_0 627 else: 628 self._eps = EPS # delattr(self, '_eps') 629 self._cD = PI_4 630 self._cE = self._cK = PI_2 631 self._cKE = _0_0 # k2 * self._cD 632 633 if alpha2: 634 if alphap2: 635 # <https://DLMF.NIST.gov/19.25.E2> 636 if kp2: 637 rj_3 = _RJ_3(_0_0, kp2, _1_0, alphap2) 638 # G(alpha2, k) 639 self._cG = self.cK + (alpha2 - k2) * rj_3 640 # H(alpha2, k) 641 self._cH = self.cK - alphap2 * rj_3 642 # Pi(alpha2, k) 643 self._cPi = self.cK + alpha2 * rj_3 644 else: 645 self._cG = self._cH = _RC(_1_0, alphap2) 646 self._cPi = INF # XXX or NAN? 647 else: 648 self._cG = self._cH = self._cPi = INF # XXX or NAN? 649 else: 650 self._cG = self.cE 651 self._cPi = self.cK 652 # cH = cK - cD but this involves large cancellations 653 # if k2 is close to 1. So write (for alpha2 = 0) 654 # cH = int(cos(phi)**2/sqrt(1-k2*sin(phi)**2),phi,0,pi/2) 655 # = 1/sqrt(1-k2) * int(sin(phi)**2/sqrt(1-k2/kp2*sin(phi)**2,...) 656 # = 1/kp * D(i*k/kp) 657 # and use D(k) = RD(0, kp2, 1) / 3 658 # so cH = 1/kp * RD(0, 1/kp2, 1) / 3 659 # = kp2 * RD(0, 1, kp2) / 3 660 # using <https://DLMF.NIST.gov/19.20.E18> 661 # Equivalently 662 # RF(x, 1) - RD(0, x, 1)/3 = x * RD(0, 1, x)/3 for x > 0 663 # For k2 = 1 and alpha2 = 0, we have 664 # cH = int(cos(phi),...) = 1 665 self._cH = kp2 * _RD_3(_0_0, _1_0, kp2) if kp2 else _1_0 666 667# self._iteration = 0 668 669 def sncndn(self, x): # PYCHOK x used! 670 '''The Jacobi elliptic function. 671 672 @arg x: The argument (C{float}). 673 674 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with 675 C{*n(B{x}, k)}. 676 677 @raise EllipticError: No convergence. 678 ''' 679 # Bulirsch's sncndn routine, p 89. 680 mc = self.kp2 681 if mc: # never negative ... 682 if mc < 0: # PYCHOK no cover 683 d = _1_0 - mc 684 mc = neg(mc / d) # /= d chokes PyChecker 685 d = sqrt(d) 686 x *= d 687 else: 688 d = 0 689 a, mn = _1_0, [] 690 for self._iteration in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC 691 # This converges quadratically, max 6 trips 692 mc = sqrt(mc) # XXX sqrt0? 693 mn.append((a, mc)) 694 c = (a + mc) * _0_5 695 if abs(a - mc) <= (_TolJAC * a): 696 break 697 mc *= a 698 a = c 699 else: 700 raise _convergenceError(self.sncndn, x) 701 x *= c 702 dn = _1_0 703 sn, cn = sincos2(x) 704 if sn: 705 a = cn / sn 706 c *= a 707 while mn: 708 m, n = mn.pop() 709 a *= c 710 c *= dn 711 dn = (n + a) / (m + a) 712 a = c / m 713 sn = copysign0(_1_0 / hypot1(c), sn) 714 cn = c * sn 715 if d: # PYCHOK no cover 716 cn, dn = dn, cn 717 sn = sn / d # /= d chokes PyChecker 718 else: 719 self._iteration = 0 720 sn = tanh(x) 721 cn = dn = _1_0 / cosh(x) 722 723 r = Elliptic3Tuple(sn, cn, dn) 724 r._iteration = self._iteration 725 return r 726 727 def _sncndn3(self, phi): 728 '''(INTERNAL) Helper for C{.fEinv} and C{._fXf}. 729 ''' 730 sn, cn = sincos2(phi) 731 return Elliptic3Tuple(sn, cn, self.fDelta(sn, cn)) 732 733 734def _convergenceError(where, *args): # PYCHOK no cover 735 '''(INTERNAL) Return an L{EllipticError}. 736 ''' 737 t = _SPACE_(where.__name__, repr(args)) 738 return EllipticError(_no_(_convergence_), txt=t) # unstr 739 740 741def _horner(e0, e1, e2, e3, e4, e5): 742 '''(INTERNAL) Horner form for C{_RD} and C{_RJ} below. 743 ''' 744 # Polynomial is <https://DLMF.NIST.gov/19.36.E2> 745 # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26 746 # - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20 + 45*E2**2*E3/272 747 # - 9*(E3*E4+E2*E5)/68) 748 # converted to Horner form ... 749 H = Fsum(471240, -540540 * e2) * e5 750 H += Fsum(612612 * e2, -540540 * e3, -556920) * e4 751 H += Fsum(306306 * e3, 675675 * e2**2, -706860 * e2, 680680) * e3 752 H += Fsum(417690 * e2, -255255 * e2**2, -875160) * e2 753 e = 4084080 * e1 754 return H.fsum_(4084080, e * e0) / e 755 756 757def _invokationError(name, *args): # PYCHOK no cover 758 '''(INTERNAL) Return an L{EllipticError}. 759 ''' 760 return EllipticError(_SPACE_('invokation', NN(name, repr(args)))) # unstr 761 762 763def _Q(A, T, tol): 764 '''(INTERNAL) Helper for C{_RD}, C{_RF} and C{_RJ}. 765 ''' 766 return max(abs(A - t) for t in T[1:]) / tol 767 768 769def _RC(x, y): # used by testElliptic.py 770 '''Degenerate symmetric integral of the first kind C{_RC(x, y)}. 771 772 @return: C{_RC(x, y) = _RF(x, y, y)}. 773 774 @see: U{C{_RC} definition<https://DLMF.NIST.gov/19.2.E17>}. 775 ''' 776 # Defined only for y != 0 and x >= 0. 777 d = x - y 778 if d < 0: # catch _NaN 779 # <https://DLMF.NIST.gov/19.2.E18> 780 d = -d 781 r = atan(sqrt(d / x)) if x > 0 else PI_2 782 elif d == _0_0: # XXX d < EPS0? 783 r, d = _1_0, y 784 elif y > 0: # <https://DLMF.NIST.gov/19.2.E19> 785 r = asinh(sqrt(d / y)) # atanh(sqrt((x - y) / x)) 786 elif y < 0: # <https://DLMF.NIST.gov/19.2.E20> 787 r = asinh(sqrt(-x / y)) # atanh(sqrt(x / (x - y))) 788 else: 789 raise _invokationError(_RC.__name__, x, y) 790 return r / sqrt(d) 791 792 793def _RD_3(x, y, z): 794 '''Degenerate symmetric integral of the third kind C{_RD(x, y, z) / 3}. 795 ''' 796 return _RD(x, y, z) / _3_0 797 798 799def _RD(x, y, z): # used by testElliptic.py 800 '''Degenerate symmetric integral of the third kind C{_RD(x, y, z)}. 801 802 @return: C{_RD(x, y, z) = _RJ(x, y, z, z)}. 803 804 @see: U{C{_RD} definition<https://DLMF.NIST.gov/19.16.E5>} and 805 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 806 ''' 807 # Carlson, eqs 2.28 - 2.34 808 A = fsum_(x, y, _3_0 * z) / _5_0 809 T = (A, x, y, z) 810 Q = _Q(A, T, _TolRD) 811 S = Fsum() 812 m = _1_0 813 for _ in range(_TRIPS): 814 An = T[0] 815 if Q < abs(m * An): # max 7 trips 816 break 817 t = T[3] # z0 818 r, s, T = _rsT3(T) 819 S += _1_0 / (m * s[2] * (t + r)) 820 m *= _4_0 821 else: 822 raise _convergenceError(_RD, x, y, z) 823 824 m *= An 825 x = (x - A) / m 826 y = (y - A) / m 827 z = (x + y) / _3_0 828 z2 = z**2 829 xy = x * y 830 return _horner(3 * S.fsum(), m * sqrt(An), 831 xy - _6_0 * z2, 832 (xy * _3_0 - _8_0 * z2) * z, 833 (xy - z2) * _3_0 * z2, 834 xy * z2 * z) 835 836 837def _RF_(x, y): 838 '''Symmetric integral of the first kind C{_RF_(x, y)}. 839 840 @return: C{_RF(x, y)}. 841 842 @see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>}. 843 ''' 844 # Carlson, eqs 2.36 - 2.38 845 a, b = sqrt(x), sqrt(y) 846 if a < b: 847 a, b = b, a 848 for _ in range(_TRIPS): 849 if abs(a - b) <= (_TolRG0 * a): # max 4 trips 850 return PI / (a + b) 851 b, a = sqrt(a * b), (a + b) * _0_5 852 853 raise _convergenceError(_RF_, x, y) 854 855 856def _RF(x, y, z): # used by testElliptic.py 857 '''Symmetric integral of the first kind C{_RF(x, y, z)}. 858 859 @return: C{_RF(x, y, z)}. 860 861 @see: U{C{_RF} definition<https://DLMF.NIST.gov/19.16.E1>} and 862 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 863 ''' 864 # Carlson, eqs 2.2 - 2.7 865 A = fmean_(x, y, z) 866 T = (A, x, y, z) 867 Q = _Q(A, T, _TolRF) 868 m = _1_0 869 for _ in range(_TRIPS): 870 An = T[0] 871 if Q < abs(m * An): # max 6 trips 872 break 873 _, _, T = _rsT3(T) 874 m *= _4_0 875 else: 876 raise _convergenceError(_RF, x, y, z) 877 878 m *= An 879 x = (A - x) / m 880 y = (A - y) / m 881 z = neg(x + y) 882 883 e2 = x * y - z**2 884 e3 = x * y * z 885 # Polynomial is <https://DLMF.NIST.gov/19.36.E1> 886 # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44 887 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16) 888 # converted to Horner form ... 889 H = Fsum( 6930 * e3, 15015 * e2**2, -16380 * e2, 17160) * e3 890 H += Fsum(10010 * e2, -5775 * e2**2, -24024) * e2 891 return H.fsum_(240240) / (240240 * sqrt(An)) 892 893 894def _RG_(x, y): 895 '''Symmetric integral of the second kind C{_RG_(x, y)}. 896 897 @return: C{_RG(x, y)}. 898 899 @see: U{C{_RG} definition<https://DLMF.NIST.gov/19.16.E3>} and 900 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 901 ''' 902 # Carlson, eqs 2.36 - 2.39 903 a, b = sqrt(x), sqrt(y) 904 if a < b: 905 a, b = b, a 906 S = Fsum(_0_25 * (a + b)**2) 907 m = _0_5 908 for _ in range(_TRIPS): # max 4 trips 909 if abs(a - b) <= (_TolRG0 * a): 910 S *= PI_2 / (a + b) 911 return S.fsum() 912 b, a = sqrt(a * b), (a + b) * _0_5 913 S -= m * (a - b)**2 914 m *= _2_0 915 916 raise _convergenceError(_RG_, x, y) 917 918 919def _RG(x, y, z): # used by testElliptic.py 920 '''Symmetric integral of the second kind C{_RG(x, y, z)}. 921 922 @return: C{_RG(x, y, z)}. 923 924 @see: U{C{_RG} definition<https://DLMF.NIST.gov/19.16.E3>} and 925 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 926 ''' 927 if not z: 928 y, z = z, y 929 # Carlson, eq 1.7 930 return fsum_(_RF(x, y, z) * z, 931 _RD_3(x, y, z) * (x - z) * (z - y), # - (y - z) 932 sqrt(x * y / z)) * _0_5 933 934 935def _RJ_3(x, y, z, p): 936 '''Symmetric integral of the third kind C{_RJ(x, y, z, p) / 3}. 937 ''' 938 return _RJ(x, y, z, p) / _3_0 939 940 941def _RJ(x, y, z, p): # used by testElliptic.py 942 '''Symmetric integral of the third kind C{_RJ(x, y, z, p)}. 943 944 @return: C{_RJ(x, y, z, p)}. 945 946 @see: U{C{_RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and 947 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 948 ''' 949 def _xyzp(x, y, z, p): 950 return (x + p) * (y + p) * (z + p) 951 952 # Carlson, eqs 2.17 - 2.25 953 D = neg(_xyzp(x, y, z, -p)) 954 A = fsum_(x, y, z, _2_0 * p) / _5_0 955 T = (A, x, y, z, p) 956 Q = _Q(A, T, _TolRD) 957 S = Fsum() 958 m = m3 = _1_0 959 for _ in range(_TRIPS): 960 An = T[0] 961 if Q < abs(m * An): # max 7 trips 962 break 963 _, s, T = _rsT3(T) 964 d = _xyzp(*s) 965 e = D / (m3 * d**2) 966 S += _RC(_1_0, _1_0 + e) / (m * d) 967 m *= _4_0 968 m3 *= 64 969 else: 970 raise _convergenceError(_RJ, x, y, z, p) 971 972 m *= An 973 x = (A - x) / m 974 y = (A - y) / m 975 z = (A - z) / m 976 xyz = x * y * z 977 p = neg(x + y + z) * _0_5 978 p2 = p**2 979 p3 = p * p2 980 981 e2 = fsum_(x * y, x * z, y * z, -p2 * _3_0) 982 return _horner(_6_0 * S.fsum(), m * sqrt(An), e2, 983 fsum_(xyz, _2_0 * p * e2, _4_0 * p3), 984 fsum_(xyz * _2_0, p * e2, _3_0 * p3) * p, 985 p2 * xyz) 986 987 988def _rsT3(T): 989 '''(INTERNAL) Helper for C{_RD}, C{_RF} and C{_RJ}. 990 ''' 991 s = map2(sqrt, T[1:]) 992 r = fdot(s[:3], s[1], s[2], s[0]) 993 T = tuple((t + r) * _0_25 for t in T) 994 return r, s, T 995 996# **) MIT License 997# 998# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. 999# 1000# Permission is hereby granted, free of charge, to any person obtaining a 1001# copy of this software and associated documentation files (the "Software"), 1002# to deal in the Software without restriction, including without limitation 1003# the rights to use, copy, modify, merge, publish, distribute, sublicense, 1004# and/or sell copies of the Software, and to permit persons to whom the 1005# Software is furnished to do so, subject to the following conditions: 1006# 1007# The above copyright notice and this permission notice shall be included 1008# in all copies or substantial portions of the Software. 1009# 1010# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 1011# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 1012# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 1013# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 1014# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 1015# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 1016# OTHER DEALINGS IN THE SOFTWARE. 1017