Coverage for pygeodesy/elliptic.py: 99%

455 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-20 14:00 -0400

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

10 

11Python 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 

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

26 

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. 

30 

31B{Elliptic integrals and functions.} 

32 

33This provides the elliptic functions and integrals needed for 

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

35categories 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 

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

49 

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

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

52 

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

54 

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

56period I{C{pi}}. 

57 

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

64 

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

69 

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 

76 

77from pygeodesy.basics import copysign0, map2, neg 

78from pygeodesy.constants import EPS, _EPStol as _TolJAC, INF, PI, PI_2, PI_4, \ 

79 _0_0, _0_125, _0_25, _0_5, _1_0, _1_64th, _2_0, \ 

80 _N_2_0, _3_0, _4_0, _6_0, _8_0, _180_0, _360_0 

81from pygeodesy.errors import _ValueError, _xkwds_pop 

82from pygeodesy.fmath import fdot, hypot1, Fsum 

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

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

85from pygeodesy.karney import _ALL_LAZY, _signBit 

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

87from pygeodesy.named import _Named, _NamedTuple, _NotImplemented 

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

89 _update_all 

90from pygeodesy.streprs import Fmt, unstr 

91from pygeodesy.units import Scalar, Scalar_ 

92from pygeodesy.utily import sincos2, sincos2d 

93 

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

95 

96__all__ = _ALL_LAZY.elliptic 

97__version__ = '23.08.19' 

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 

106class _Complete(object): 

107 '''(INTERAL) Hold complete integrals. 

108 ''' 

109 def __init__(self, **kwds): 

110 self.__dict__ = kwds 

111 

112 

113class _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 

154class _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 

167class 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, *over): 

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

830 

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

832 / over} with C{over} typically 3. 

833 

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

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

836 ''' 

837 return _RD(None, x, y, z, *over) 

838 

839 @staticmethod 

840 def fRF(x, y, z=0): 

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

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

843 

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

845 

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

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

848 ''' 

849 return _RF3(None, x, y, z) if z else _RF2(None, x, y) 

850 

851 @staticmethod 

852 def fRG(x, y, z=0): 

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

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

855 

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

857 

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

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

860 ''' 

861 return _RG3(None, x, y, z) if z else (_RG2(None, x, y) * _0_5) 

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 

877class EllipticError(_ValueError): 

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

879 ''' 

880 pass 

881 

882 

883class 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 

890class _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 try: 

907 s = Fsum(*t).fover(n) 

908 except ValueError: # NAN 

909 s = sum(t) / n 

910 self._a0 = self._a = s 

911 return s 

912 

913 def asr3(self, a): 

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

915 ''' 

916 L = self 

917 # assert a is L._a 

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

919 try: 

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

921 except ValueError: # NAN 

922 r = sum(s[i] * s[(i + 1) % 3] for i in range(3)) 

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

924 # assert L is self 

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

926 return a, s, r 

927 

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

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

930 ''' 

931 for x in xy_: 

932 yield (self._a0 - x) / am 

933 

934 def thresh(self, Tol): 

935 '''Return the convergence threshold. 

936 ''' 

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

938 

939 

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

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

942 ''' 

943 E22 = E2**2 

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

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

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

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

948 # converted to Horner-like form ... 

949 e = e1 * 4084080 

950 S *= e 

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

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

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

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

955 S += 4084080 

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

957 

958 

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

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

961 ''' 

962 n = _DOT_(Elliptic.__name__, name) 

963 n = _SPACE_(_invokation_, n) 

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

965 

966 

967def _iterations(inst, i): 

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

969 ''' 

970 if inst: 

971 inst._iteration += i 

972 

973 

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

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

976 ''' 

977 n = Elliptic.__name__ 

978 if where: 

979 n = _DOT_(n, where.__name__) 

980 if kwds_thresh: 

981 q = _xkwds_pop(kwds_thresh, thresh=False) 

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

983 else: 

984 q = False 

985 t = unstr(n, *args) 

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

987 

988 

989def _3over(a, b): 

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

991 ''' 

992 return _3_0 / (a * b) 

993 

994 

995def _RC(unused, x, y): 

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

997 ''' 

998 d = x - y 

999 if d < 0: # catch _NaN 

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

1001 d = -d 

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

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

1004 d, r = y, _1_0 

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

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

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

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

1009 else: 

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

1011 return r / sqrt(d) 

1012 

1013 

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

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

1016 ''' 

1017 L = _Lxyz(x, y, z) 

1018 a = L.a0(5) 

1019 q = L.thresh(_TolRF) 

1020 S = _Deferred() 

1021 am, m = a, 1 

1022 for i in range(_TRIPS): 

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

1024 _iterations(inst, i) 

1025 break 

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

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

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

1029 m *= 4 

1030 am = a * m 

1031 else: # PYCHOK no cover 

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

1033 thresh=True) 

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

1035 xy = x * y 

1036 z = (x + y) / _3_0 

1037 z2 = z**2 

1038 S = S.Fsum.fmul(_3_0) 

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

1040 xy - _6_0 * z2, 

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

1042 (xy - z2) * _3_0 * z2, 

1043 xy * z2 * z, *over) 

1044 

1045 

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

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

1048 ''' 

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

1050 if a < b: 

1051 a, b = b, a 

1052 for i in range(_TRIPS): 

1053 t = _TolRG0 * a 

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

1055 _iterations(inst, i) 

1056 return (PI / (a + b)) 

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

1058 else: # PYCHOK no cover 

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

1060 

1061 

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

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

1064 ''' 

1065 L = _Lxyz(x, y, z) 

1066 a = L.a0(3) 

1067 q = L.thresh(_TolRF) 

1068 am, m = a, 1 

1069 for i in range(_TRIPS): 

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

1071 _iterations(inst, i) 

1072 break 

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

1074 m *= 4 

1075 am = a * m 

1076 else: # PYCHOK no cover 

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

1078 thresh=True) 

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

1080 z = neg(x + y) 

1081 xy = x * y 

1082 e2 = xy - z**2 

1083 e3 = xy * z 

1084 e4 = e2**2 

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

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

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

1088 # converted to Horner-like form ... 

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

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

1091 S += Fsum(240240) 

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

1093 

1094 

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

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

1097 ''' 

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

1099 if a < b: 

1100 a, b = b, a 

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

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

1103 m = -1 

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

1105 t = _TolRG0 * a 

1106 if ab <= t: 

1107 _iterations(inst, i) 

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

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

1110 ab = fabs(a - b) 

1111 S += ab**2 * m 

1112 m *= 2 

1113 else: # PYCHOK no cover 

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

1115 

1116 

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

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

1119 ''' 

1120# if not z: 

1121# y, z = z, y 

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

1123 if rd: # Carlson, eq 1.7 

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

1125 xyz = x * y 

1126 if xyz: 

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

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

1129 

1130 

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

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

1133 ''' 

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

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

1136 

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

1138 a = L.a0(5) 

1139 q = L.thresh(_TolRD) 

1140 S = _Deferred() 

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

1142 am, m = a, 1 

1143 for i in range(_TRIPS): 

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

1145 _iterations(inst, i) 

1146 break 

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

1148 d = _xyzp(*s) 

1149 if n: 

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

1151 n *= _1_64th # /= chokes PyChecker 

1152 else: 

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

1154 S += rc / (d * m) 

1155 m *= 4 

1156 am = a * m 

1157 else: # PYCHOK no cover 

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

1159 thresh=True) 

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

1161 xyz = x * y * z 

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

1163 p2 = p**2 

1164 p3 = p2*p 

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

1166 E2p = E2 * p 

1167 S = S.Fsum.fmul(_6_0) 

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

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

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

1171 xyz * p2, *over) 

1172 

1173# **) MIT License 

1174# 

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

1176# 

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

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

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

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

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

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

1183# 

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

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

1186# 

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

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

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

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

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

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

1193# OTHER DEALINGS IN THE SOFTWARE.