Coverage for pygeodesy/elliptic.py: 96%

485 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-01 11:43 -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, neg_ 

78from pygeodesy.constants import EPS, INF, NAN, PI, PI_2, PI_4, \ 

79 _EPStol as _TolJAC, _0_0, _1_64th, \ 

80 _0_25, _0_5, _1_0, _2_0, _N_2_0, \ 

81 _3_0, _4_0, _6_0, _8_0, _180_0, \ 

82 _360_0, _over 

83# from pygeodesy.errors import _ValueError # from .fmath 

84from pygeodesy.fmath import fdot, hypot1, zqrt, _ValueError 

85from pygeodesy.fsums import Fsum, _sum 

86# from pygeodesy.internals import _dunder_nameof # from .lazily 

87from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _invalid_, \ 

88 _invokation_, _negative_, _SPACE_ 

89from pygeodesy.karney import _K_2_0, _norm180, _signBit, _sincos2 

90from pygeodesy.lazily import _ALL_LAZY, _dunder_nameof 

91from pygeodesy.named import _Named, _NamedTuple, Fmt, unstr 

92from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all 

93# from pygeodesy.streprs import Fmt, unstr # from .named 

94from pygeodesy.units import Scalar, Scalar_ 

95# from pygeodesy.utily import sincos2 as _sincos2 # from .karney 

96 

97from math import asinh, atan, atan2, ceil, cosh, fabs, floor, \ 

98 radians, sin, sqrt, tanh 

99 

100__all__ = _ALL_LAZY.elliptic 

101__version__ = '24.05.18' 

102 

103_TolRD = zqrt(EPS * 0.002) 

104_TolRF = zqrt(EPS * 0.030) 

105_TolRG0 = _TolJAC * 2.7 

106_TRIPS = 21 # Max depth, 7 might be sufficient 

107 

108 

109class _Cs(object): 

110 '''(INTERAL) Complete integrals cache. 

111 ''' 

112 def __init__(self, **kwds): 

113 self.__dict__ = kwds 

114 

115 

116class _Dsum(list): 

117 '''(INTERNAL) Deferred C{Fsum}. 

118 ''' 

119 def __call__(self, s): 

120 try: # Fsum *= s 

121 return Fsum(*self).fmul(s) 

122 except ValueError: # Fsum(NAN) exception 

123 return _sum(self) * s 

124 

125 def __iadd__(self, x): 

126 list.append(self, x) 

127 return self 

128 

129 

130class Elliptic(_Named): 

131 '''Elliptic integrals and functions. 

132 

133 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/ 

134 C++/doc/classGeographicLib_1_1EllipticFunction.html#details>}. 

135 ''' 

136# _alpha2 = 0 

137# _alphap2 = 0 

138# _eps = EPS 

139# _k2 = 0 

140# _kp2 = 0 

141 

142 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None, **name): 

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

144 

145 @kwarg name: Optional C{B{name}=NN} (C{str}). 

146 

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

148 

149 @note: If only elliptic integrals of the first and second kinds 

150 are needed, use C{B{alpha2}=0}, the default value. In 

151 that case, we have C{Π(φ, 0, k) = F(φ, k), G(φ, 0, k) = 

152 E(φ, k)} and C{H(φ, 0, k) = F(φ, k) - D(φ, k)}. 

153 ''' 

154 self.reset(k2=k2, alpha2=alpha2, kp2=kp2, alphap2=alphap2) 

155 if name: 

156 self.name = name 

157 

158 @Property_RO 

159 def alpha2(self): 

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

161 ''' 

162 return self._alpha2 

163 

164 @Property_RO 

165 def alphap2(self): 

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

167 ''' 

168 return self._alphap2 

169 

170 @Property_RO 

171 def cD(self): 

172 '''Get Jahnke's complete integral C{D(k)} (C{float}), 

173 U{defined<https://DLMF.NIST.gov/19.2.E6>}. 

174 ''' 

175 return self._cDEKEeps.cD 

176 

177 @Property_RO 

178 def _cDEKEeps(self): 

179 '''(INTERNAL) Get the complete integrals D, E, K and KE plus C{eps}. 

180 ''' 

181 k2, kp2 = self.k2, self.kp2 

182 if k2: 

183 if kp2: 

184 try: 

185 self._iteration = 0 

186 # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3 

187 # <https://DLMF.NIST.gov/19.25.E1> 

188 D = _RD(self, _0_0, kp2, _1_0, _3_0) 

189 cD = float(D) 

190 # Complete elliptic integral E(k), Carlson eq. 4.2 

191 # <https://DLMF.NIST.gov/19.25.E1> 

192 cE = _rG2(self, kp2, _1_0, PI_=PI_2) 

193 # Complete elliptic integral K(k), Carlson eq. 4.1 

194 # <https://DLMF.NIST.gov/19.25.E1> 

195 cK = _rF2(self, kp2, _1_0) 

196 cKE = float(D.fmul(k2)) 

197 eps = k2 / (sqrt(kp2) + _1_0)**2 

198 

199 except Exception as e: 

200 raise _ellipticError(self.reset, k2=k2, kp2=kp2, cause=e) 

201 else: 

202 cD = cK = cKE = INF 

203 cE = _1_0 

204 eps = k2 

205 else: 

206 cD = PI_4 

207 cE = cK = PI_2 

208 cKE = _0_0 # k2 * cD 

209 eps = EPS 

210 

211 return _Cs(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps) 

212 

213 @Property_RO 

214 def cE(self): 

215 '''Get the complete integral of the second kind C{E(k)} 

216 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}. 

217 ''' 

218 return self._cDEKEeps.cE 

219 

220 @Property_RO 

221 def cG(self): 

222 '''Get Legendre's complete geodesic longitude integral 

223 C{G(α^2, k)} (C{float}). 

224 ''' 

225 return self._cGHPi.cG 

226 

227 @Property_RO 

228 def _cGHPi(self): 

229 '''(INTERNAL) Get the complete integrals G, H and Pi. 

230 ''' 

231 alpha2, alphap2, kp2 = self.alpha2, self.alphap2, self.kp2 

232 try: 

233 self._iteration = 0 

234 if alpha2: 

235 if alphap2: 

236 if kp2: # <https://DLMF.NIST.gov/19.25.E2> 

237 cK = self.cK 

238 Rj = _RJ(self, _0_0, kp2, _1_0, alphap2, _3_0) 

239 cG = float(Rj * (alpha2 - self.k2) + cK) # G(alpha2, k) 

240 cH = -float(Rj * alphap2 - cK) # H(alpha2, k) 

241 cPi = float(Rj * alpha2 + cK) # Pi(alpha2, k) 

242 else: 

243 cG = cH = _rC(self, _1_0, alphap2) 

244 cPi = INF # XXX or NAN? 

245 else: 

246 cG = cH = cPi = INF # XXX or NAN? 

247 else: 

248 cG, cPi = self.cE, self.cK 

249 # H = K - D but this involves large cancellations if k2 is near 1. 

250 # So write (for alpha2 = 0) 

251 # H = int(cos(phi)**2 / sqrt(1-k2 * sin(phi)**2), phi, 0, pi/2) 

252 # = 1 / sqrt(1-k2) * int(sin(phi)**2 / sqrt(1-k2/kp2 * sin(phi)**2,...) 

253 # = 1 / kp * D(i * k/kp) 

254 # and use D(k) = RD(0, kp2, 1) / 3, so 

255 # H = 1/kp * RD(0, 1/kp2, 1) / 3 

256 # = kp2 * RD(0, 1, kp2) / 3 

257 # using <https://DLMF.NIST.gov/19.20.E18>. Equivalently 

258 # RF(x, 1) - RD(0, x, 1) / 3 = x * RD(0, 1, x) / 3 for x > 0 

259 # For k2 = 1 and alpha2 = 0, we have 

260 # H = int(cos(phi),...) = 1 

261 cH = float(_RD(self, _0_0, _1_0, kp2, _3_0 / kp2)) if kp2 else _1_0 

262 

263 except Exception as e: 

264 raise _ellipticError(self.reset, kp2=kp2, alpha2 =alpha2, 

265 alphap2=alphap2, cause=e) 

266 return _Cs(cG=cG, cH=cH, cPi=cPi) 

267 

268 @Property_RO 

269 def cH(self): 

270 '''Get Cayley's complete geodesic longitude difference integral 

271 C{H(α^2, k)} (C{float}). 

272 ''' 

273 return self._cGHPi.cH 

274 

275 @Property_RO 

276 def cK(self): 

277 '''Get the complete integral of the first kind C{K(k)} 

278 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}. 

279 ''' 

280 return self._cDEKEeps.cK 

281 

282 @Property_RO 

283 def cKE(self): 

284 '''Get the difference between the complete integrals of the 

285 first and second kinds, C{K(k) − E(k)} (C{float}). 

286 ''' 

287 return self._cDEKEeps.cKE 

288 

289 @Property_RO 

290 def cPi(self): 

291 '''Get the complete integral of the third kind C{Pi(α^2, k)} 

292 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}. 

293 ''' 

294 return self._cGHPi.cPi 

295 

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

297 '''Jahnke's periodic incomplete elliptic integral. 

298 

299 @arg sn: sin(φ). 

300 @arg cn: cos(φ). 

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

302 

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

304 

305 @raise EllipticError: Invalid invokation or no convergence. 

306 ''' 

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

308 

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

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

311 

312 @arg sn: sin(φ). 

313 @arg cn: cos(φ). 

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

315 

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

317 

318 @raise EllipticError: Invalid invokation or no convergence. 

319 ''' 

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

321 

322 def deltaEinv(self, stau, ctau): 

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

324 

325 @arg stau: sin(τ) 

326 @arg ctau: cos(τ) 

327 

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

329 

330 @raise EllipticError: No convergence. 

331 ''' 

332 try: 

333 if _signBit(ctau): # pi periodic 

334 stau, ctau = neg_(stau, ctau) 

335 t = atan2(stau, ctau) 

336 return self._Einv(t * self.cE / PI_2) - t 

337 

338 except Exception as e: 

339 raise _ellipticError(self.deltaEinv, stau, ctau, cause=e) 

340 

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

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

343 

344 @arg sn: sin(φ). 

345 @arg cn: cos(φ). 

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

347 

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

349 

350 @raise EllipticError: Invalid invokation or no convergence. 

351 ''' 

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

353 

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

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

356 

357 @arg sn: sin(φ). 

358 @arg cn: cos(φ). 

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

360 

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

362 

363 @raise EllipticError: Invalid invokation or no convergence. 

364 ''' 

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

366 

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

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

369 

370 @arg sn: sin(φ). 

371 @arg cn: cos(φ). 

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

373 

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

375 

376 @raise EllipticError: Invalid invokation or no convergence. 

377 ''' 

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

379 

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

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

382 

383 @arg sn: sin(φ). 

384 @arg cn: cos(φ). 

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

386 

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

388 (C{float}). 

389 

390 @raise EllipticError: Invalid invokation or no convergence. 

391 ''' 

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

393 

394 def _Einv(self, x): 

395 '''(INTERNAL) Helper for C{.deltaEinv} and C{.fEinv}. 

396 ''' 

397 E2 = self.cE * _2_0 

398 n = floor(x / E2 + _0_5) 

399 r = x - E2 * n # r in [-cE, cE) 

400 # linear approximation 

401 phi = PI * r / E2 # phi in [-PI_2, PI_2) 

402 Phi = Fsum(phi) 

403 # first order correction 

404 phi = Phi.fsum_(self.eps * sin(phi * _2_0) / _N_2_0) 

405 self._iteration = 0 

406 # For kp2 close to zero use asin(r / cE) or J. P. Boyd, 

407 # Applied Math. and Computation 218, 7005-7013 (2012) 

408 # <https://DOI.org/10.1016/j.amc.2011.12.021> 

409 _Phi2 = Phi.fsum2f_ # aggregate 

410 _abs = fabs 

411 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC 

412 sn, cn, dn = self._sncndn3(phi) 

413 if dn: 

414 sn = self.fE(sn, cn, dn) 

415 phi, d = _Phi2((r - sn) / dn) 

416 else: # PYCHOK no cover 

417 d = _0_0 # XXX continue? 

418 if _abs(d) < _TolJAC: # 3-4 trips 

419 _iterations(self, i) 

420 break 

421 else: # PYCHOK no cover 

422 raise _convergenceError(d, _TolJAC) 

423 return Phi.fsum_(n * PI) if n else phi 

424 

425 @Property_RO 

426 def eps(self): 

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

428 ''' 

429 return self._cDEKEeps.eps 

430 

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

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

433 Jacobi elliptic functions. 

434 

435 @arg phi_or_sn: φ or sin(φ). 

436 @kwarg cn: C{None} or cos(φ). 

437 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

438 

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

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

441 

442 @raise EllipticError: Invalid invokation or no convergence. 

443 ''' 

444 def _fD(sn, cn, dn): 

445 r = fabs(sn)**3 

446 if r: 

447 r = float(_RD(self, cn**2, dn**2, _1_0, _3_0 / r)) 

448 return r 

449 

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

451 self.deltaD, _fD) 

452 

453 def fDelta(self, sn, cn): 

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

455 

456 @arg sn: sin(φ). 

457 @arg cn: cos(φ). 

458 

459 @return: sqrt(1 − k2 * sin(2φ)) (C{float}). 

460 ''' 

461 try: 

462 k2 = self.k2 

463 s = (self.kp2 + cn**2 * k2) if k2 > 0 else ( 

464 (_1_0 - sn**2 * k2) if k2 < 0 else self.kp2) 

465 return sqrt(s) if s else _0_0 

466 

467 except Exception as e: 

468 raise _ellipticError(self.fDelta, sn, cn, k2=k2, cause=e) 

469 

470 def fE(self, phi_or_sn, cn=None, dn=None): 

471 '''The incomplete integral of the second kind in terms of 

472 Jacobi elliptic functions. 

473 

474 @arg phi_or_sn: φ or sin(φ). 

475 @kwarg cn: C{None} or cos(φ). 

476 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

477 

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

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

480 

481 @raise EllipticError: Invalid invokation or no convergence. 

482 ''' 

483 def _fE(sn, cn, dn): 

484 '''(INTERNAL) Core of C{.fE}. 

485 ''' 

486 if sn: 

487 sn2, cn2, dn2 = sn**2, cn**2, dn**2 

488 kp2, k2 = self.kp2, self.k2 

489 if k2 <= 0: # Carlson, eq. 4.6, <https://DLMF.NIST.gov/19.25.E9> 

490 Ei = _RF3(self, cn2, dn2, _1_0) 

491 if k2: 

492 Ei -= _RD(self, cn2, dn2, _1_0, _3over(k2, sn2)) 

493 elif kp2 >= 0: # k2 > 0, <https://DLMF.NIST.gov/19.25.E10> 

494 Ei = _over(k2 * fabs(cn), dn) # float 

495 if kp2: 

496 Ei += (_RD( self, cn2, _1_0, dn2, _3over(k2, sn2)) + 

497 _RF3(self, cn2, dn2, _1_0)) * kp2 

498 else: # kp2 < 0, <https://DLMF.NIST.gov/19.25.E11> 

499 Ei = _over(dn, fabs(cn)) 

500 Ei -= _RD(self, dn2, _1_0, cn2, _3over(kp2, sn2)) 

501 Ei *= fabs(sn) 

502 ei = float(Ei) 

503 else: # PYCHOK no cover 

504 ei = _0_0 

505 return ei 

506 

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

508 self.deltaE, _fE) 

509 

510 def fEd(self, deg): 

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

512 the argument given in C{degrees}. 

513 

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

515 

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

517 

518 @raise EllipticError: No convergence. 

519 ''' 

520 if _K_2_0: 

521 e = round((deg - _norm180(deg)) / _360_0) 

522 elif fabs(deg) < _180_0: 

523 e = _0_0 

524 else: 

525 e = ceil(deg / _360_0 - _0_5) 

526 deg -= e * _360_0 

527 return self.fE(radians(deg)) + e * self.cE * _4_0 

528 

529 def fEinv(self, x): 

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

531 

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

533 

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

535 (C{float}). 

536 

537 @raise EllipticError: No convergence. 

538 ''' 

539 try: 

540 return self._Einv(x) 

541 except Exception as e: 

542 raise _ellipticError(self.fEinv, x, cause=e) 

543 

544 def fF(self, phi_or_sn, cn=None, dn=None): 

545 '''The incomplete integral of the first kind in terms of 

546 Jacobi elliptic functions. 

547 

548 @arg phi_or_sn: φ or sin(φ). 

549 @kwarg cn: C{None} or cos(φ). 

550 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

551 

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

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

554 

555 @raise EllipticError: Invalid invokation or no convergence. 

556 ''' 

557 def _fF(sn, cn, dn): 

558 r = fabs(sn) 

559 if r: 

560 r = float(_RF3(self, cn**2, dn**2, _1_0).fmul(r)) 

561 return r 

562 

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

564 self.deltaF, _fF) 

565 

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

567 '''Legendre's geodesic longitude integral 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: G(φ, k) as though φ ∈ (−π, π] (C{float}). 

575 

576 @raise EllipticError: Invalid invokation or no convergence. 

577 

578 @note: Legendre expresses the longitude of a point on the 

579 geodesic in terms of this combination of elliptic 

580 integrals in U{Exercices de Calcul Intégral, Vol 1 

581 (1811), p 181<https://Books.Google.com/books?id= 

582 riIOAAAAQAAJ&pg=PA181>}. 

583 

584 @see: U{Geodesics in terms of elliptic integrals<https:// 

585 GeographicLib.SourceForge.io/C++/doc/geodesic.html#geodellip>} 

586 for the expression for the longitude in terms of this function. 

587 ''' 

588 return self._fXa(phi_or_sn, cn, dn, self.alpha2 - self.k2, 

589 self.cG, self.deltaG) 

590 

591 def fH(self, phi_or_sn, cn=None, dn=None): 

592 '''Cayley's geodesic longitude difference integral in terms of 

593 Jacobi elliptic functions. 

594 

595 @arg phi_or_sn: φ or sin(φ). 

596 @kwarg cn: C{None} or cos(φ). 

597 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

598 

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

600 

601 @raise EllipticError: Invalid invokation or no convergence. 

602 

603 @note: Cayley expresses the longitude difference of a point 

604 on the geodesic in terms of this combination of 

605 elliptic integrals in U{Phil. Mag. B{40} (1870), p 333 

606 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}. 

607 

608 @see: U{Geodesics in terms of elliptic integrals<https:// 

609 GeographicLib.SourceForge.io/C++/doc/geodesic.html#geodellip>} 

610 for the expression for the longitude in terms of this function. 

611 ''' 

612 return self._fXa(phi_or_sn, cn, dn, -self.alphap2, 

613 self.cH, self.deltaH) 

614 

615 def fPi(self, phi_or_sn, cn=None, dn=None): 

616 '''The incomplete integral of the third kind in terms of 

617 Jacobi elliptic functions. 

618 

619 @arg phi_or_sn: φ or sin(φ). 

620 @kwarg cn: C{None} or cos(φ). 

621 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

622 

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

624 

625 @raise EllipticError: Invalid invokation or no convergence. 

626 ''' 

627 if dn is None and cn is not None: # and isscalar(phi_or_sn) 

628 dn = self.fDelta(phi_or_sn, cn) # in .triaxial 

629 return self._fXa(phi_or_sn, cn, dn, self.alpha2, 

630 self.cPi, self.deltaPi) 

631 

632 def _fXa(self, phi_or_sn, cn, dn, aX, cX, deltaX): 

633 '''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}. 

634 ''' 

635 def _fX(sn, cn, dn): 

636 if sn: 

637 cn2, dn2 = cn**2, dn**2 

638 R = _RF3(self, cn2, dn2, _1_0) 

639 if aX: 

640 sn2 = sn**2 

641 p = sn2 * self.alphap2 + cn2 

642 R += _RJ(self, cn2, dn2, _1_0, p, _3over(aX, sn2)) 

643 R *= fabs(sn) 

644 r = float(R) 

645 else: # PYCHOK no cover 

646 r = _0_0 

647 return r 

648 

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

650 

651 def _fXf(self, phi_or_sn, cn, dn, cX, deltaX, fX): 

652 '''(INTERNAL) Helper for C{.fD}, C{.fE}, C{.fF} and C{._fXa}. 

653 ''' 

654 self._iteration = 0 # aggregate 

655 phi = sn = phi_or_sn 

656 if cn is dn is None: # fX(phi) call 

657 sn, cn, dn = self._sncndn3(phi) 

658 if fabs(phi) >= PI: 

659 return (deltaX(sn, cn, dn) + phi) * cX / PI_2 

660 # fall through 

661 elif cn is None or dn is None: 

662 n = NN(_f_, deltaX.__name__[5:]) 

663 raise _ellipticError(n, sn, cn, dn) 

664 

665 if _signBit(cn): # enforce usual trig-like symmetries 

666 xi = cX * _2_0 - fX(sn, cn, dn) 

667 else: 

668 xi = fX(sn, cn, dn) if cn > 0 else cX 

669 return copysign0(xi, sn) 

670 

671 @Property_RO 

672 def k2(self): 

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

674 ''' 

675 return self._k2 

676 

677 @Property_RO 

678 def kp2(self): 

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

680 ''' 

681 return self._kp2 

682 

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

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

685 

686 @kwarg k2: Modulus squared (C{float}, NINF <= k^2 <= 1). 

687 @kwarg alpha2: Parameter squared (C{float}, NINF <= α^2 <= 1). 

688 @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0). 

689 @kwarg alphap2: Complementary parameter squared (C{float}, α'^2 >= 0). 

690 

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

692 or B{C{alphap2}}. 

693 

694 @note: The arguments must satisfy C{B{k2} + B{kp2} = 1} and 

695 C{B{alpha2} + B{alphap2} = 1}. No checking is done 

696 that these conditions are met to enable accuracy to be 

697 maintained, e.g., when C{k} is very close to unity. 

698 ''' 

699 if self.__dict__: 

700 _update_all(self, _Named.iteration._uname, Base=Property_RO) 

701 

702 self._k2 = Scalar_(k2=k2, Error=EllipticError, low=None, high=_1_0) 

703 self._kp2 = Scalar_(kp2=((_1_0 - k2) if kp2 is None else kp2), Error=EllipticError) 

704 

705 self._alpha2 = Scalar_(alpha2=alpha2, Error=EllipticError, low=None, high=_1_0) 

706 self._alphap2 = Scalar_(alphap2=((_1_0 - alpha2) if alphap2 is None else alphap2), 

707 Error=EllipticError) 

708 

709 # Values of complete elliptic integrals for k = 0,1 and alpha = 0,1 

710 # K E D 

711 # k = 0: pi/2 pi/2 pi/4 

712 # k = 1: inf 1 inf 

713 # Pi G H 

714 # k = 0, alpha = 0: pi/2 pi/2 pi/4 

715 # k = 1, alpha = 0: inf 1 1 

716 # k = 0, alpha = 1: inf inf pi/2 

717 # k = 1, alpha = 1: inf inf inf 

718 # 

719 # G(0, k) = Pi(0, k) = H(1, k) = E(k) 

720 # H(0, k) = K(k) - D(k) 

721 # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2)) 

722 # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1)) 

723 # Pi(alpha2, 1) = inf 

724 # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2) 

725 

726 def sncndn(self, x): 

727 '''The Jacobi elliptic function. 

728 

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

730 

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

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

733 

734 @raise EllipticError: No convergence. 

735 ''' 

736 self._iteration = 0 # reset 

737 try: # Bulirsch's sncndn routine, p 89. 

738 if self.kp2: 

739 c, d, cd, mn = self._sncndn4 

740 dn = _1_0 

741 sn, cn = _sincos2(x * cd) 

742 if sn: 

743 a = cn / sn 

744 c *= a 

745 for m, n in reversed(mn): 

746 a *= c 

747 c *= dn 

748 dn = (n + a) / (m + a) 

749 a = c / m 

750 a = _1_0 / hypot1(c) 

751 sn = neg(a) if _signBit(sn) else a 

752 cn = c * sn 

753 if d and _signBit(self.kp2): 

754 cn, dn = dn, cn 

755 sn = sn / d # /= chokes PyChecker 

756 else: 

757 sn = tanh(x) 

758 cn = dn = _1_0 / cosh(x) 

759 

760 except Exception as e: 

761 raise _ellipticError(self.sncndn, x, kp2=self.kp2, cause=e) 

762 

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

764 

765 def _sncndn3(self, phi): 

766 '''(INTERNAL) Helper for C{.fEinv} and C{._fXf}. 

767 ''' 

768 sn, cn = _sincos2(phi) 

769 return sn, cn, self.fDelta(sn, cn) 

770 

771 @Property_RO 

772 def _sncndn4(self): 

773 '''(INTERNAL) Get Bulirsch' 4-tuple C{(c, d, cd, mn)}. 

774 ''' 

775 # Bulirsch's sncndn routine, p 89. 

776 d, mc = 0, self.kp2 

777 if _signBit(mc): 

778 d = _1_0 - mc 

779 mc = neg(mc / d) 

780 d = sqrt(d) 

781 

782 mn, a = [], _1_0 

783 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC 

784 mc = sqrt(mc) 

785 mn.append((a, mc)) 

786 c = (a + mc) * _0_5 

787 r = fabs(mc - a) 

788 t = _TolJAC * a 

789 if r <= t: # 6 trips, quadratic 

790 _iterations(self, i) 

791 break 

792 mc *= a 

793 a = c 

794 else: # PYCHOK no cover 

795 raise _convergenceError(r, t) 

796 cd = (c * d) if d else c 

797 return c, d, cd, mn 

798 

799 @staticmethod 

800 def fRC(x, y): 

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

802 

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

804 

805 @see: U{C{RC} definition<https://DLMF.NIST.gov/19.2.E17>} and 

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

807 ''' 

808 return _rC(None, x, y) 

809 

810 @staticmethod 

811 def fRD(x, y, z, *over): 

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

813 

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

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

816 

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

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

819 ''' 

820 try: 

821 return float(_RD(None, x, y, z, *over)) 

822 except Exception as e: 

823 raise _ellipticError(Elliptic.fRD, x, y, z, *over, cause=e) 

824 

825 @staticmethod 

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

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

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

829 

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

831 

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

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

834 ''' 

835 try: 

836 return float(_RF3(None, x, y, z)) if z else _rF2(None, x, y) 

837 except Exception as e: 

838 raise _ellipticError(Elliptic.fRF, x, y, z, cause=e) 

839 

840 @staticmethod 

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

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

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

844 

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

846 

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

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

849 U{RG<https://GeographicLib.SourceForge.io/C++/doc/ 

850 EllipticFunction_8cpp_source.html#l00096>} version 2.3. 

851 ''' 

852 try: 

853 return _rG2(None, x, y) if z == 0 else ( 

854 _rG2(None, z, x) if y == 0 else ( 

855 _rG2(None, y, z) if x == 0 else _rG3(None, x, y, z))) 

856 except Exception as e: 

857 t = _negative_ if min(x, y, z) < 0 else NN 

858 raise _ellipticError(Elliptic.fRG, x, y, z, cause=e, txt=t) 

859 

860 @staticmethod 

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

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

863 

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

865 

866 @see: U{C{RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and 

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

868 ''' 

869 try: 

870 return float(_RJ(None, x, y, z, p)) 

871 except Exception as e: 

872 raise _ellipticError(Elliptic.fRJ, x, y, z, p, cause=e) 

873 

874 @staticmethod 

875 def _RFRD(x, y, z, m): 

876 # in .auxilats.AuxDLat.DE, .auxilats.AuxLat.Rectifying 

877 try: # float(RF(x, y, z) - RD(x, y, z, 3 / m)) 

878 R = _RF3(None, x, y, z) 

879 if m: 

880 R -= _RD(None, x, y, z, _3_0 / m) 

881 except Exception as e: 

882 raise _ellipticError(Elliptic._RFRD, x, y, z, m, cause=e) 

883 return float(R) 

884 

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

886 

887 

888class EllipticError(_ValueError): 

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

890 ''' 

891 pass 

892 

893 

894class Elliptic3Tuple(_NamedTuple): 

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

896 ''' 

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

898 _Units_ = ( Scalar, Scalar, Scalar) 

899 

900 

901class _List(list): 

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

903 ''' 

904 _a0 = None 

905# _xyzp = () 

906 

907 def __init__(self, *xyzp): # x, y, z [, p] 

908 list.__init__(self, xyzp) 

909 self._xyzp = xyzp 

910 

911 def a0(self, n): 

912 '''Compute the initial C{a}. 

913 ''' 

914 t = tuple(self) 

915 m = n - len(t) 

916 if m > 0: 

917 t += t[-1:] * m 

918 try: 

919 a = Fsum(*t).fover(n) 

920 except ValueError: # Fsum(NAN) exception 

921 a = _sum(t) / n 

922 self._a0 = a 

923 return a 

924 

925 def amrs4(self, inst, y, Tol): 

926 '''Yield Carlson 4-tuples C{(An, mul, lam, s)} plus sentinel, with 

927 C{lam = fdot(sqrt(x), ... (z))} and C{s = (sqrt(x), ... (p))}. 

928 ''' 

929 L = self 

930 a = L.a0(5 if y else 3) 

931 m = 1 

932 t = max(fabs(a - _) for _ in L) / Tol 

933 for i in range(_TRIPS): 

934 d = fabs(a * m) 

935 if d > t: # 3-6 trips 

936 _iterations(inst, i) 

937 break 

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

939 try: 

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

941 except ValueError: # Fsum(NAN) exception 

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

943 L[:] = ((r + _) * _0_25 for _ in L) 

944 a = (r + a) * _0_25 

945 if y: # yield only if used 

946 yield a, m, r, s # L[2] is next z 

947 m *= 4 

948 else: # PYCHOK no cover 

949 raise _convergenceError(d, t, thresh=True) 

950 yield a, m, None, () # sentinel: same a, next m, no r and s 

951 

952 def rescale(self, am, *xs): 

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

954 ''' 

955 # assert am 

956 a0 = self._a0 

957 _am = _1_0 / am 

958 for x in xs: 

959 yield (a0 - x) * _am 

960 

961 

962def _ab2(inst, x, y): 

963 '''(INTERNAL) Yield Carlson 2-tuples C{(xn, yn)}. 

964 ''' 

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

966 if b > a: 

967 b, a = a, b 

968 for i in range(_TRIPS): 

969 yield a, b # xi, yi 

970 d = fabs(a - b) 

971 t = _TolRG0 * a 

972 if d <= t: # 3-4 trips 

973 _iterations(inst, i) 

974 break 

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

976 else: # PYCHOK no cover 

977 raise _convergenceError(d, t) 

978 

979 

980def _convergenceError(d, tol, **thresh): 

981 '''(INTERNAL) Format a no-convergence Error. 

982 ''' 

983 t = Fmt.no_convergence(d, tol, **thresh) 

984 return ValueError(t) # txt only 

985 

986 

987def _deltaX(sn, cn, dn, cX, fX): 

988 '''(INTERNAL) Helper for C{Elliptic.deltaD} thru C{.deltaPi}. 

989 ''' 

990 try: 

991 if cn is None or dn is None: 

992 raise ValueError(_invalid_) 

993 

994 if _signBit(cn): 

995 sn, cn = neg_(sn, cn) 

996 r = fX(sn, cn, dn) * PI_2 / cX 

997 return r - atan2(sn, cn) 

998 

999 except Exception as e: 

1000 n = NN(_delta_, fX.__name__[1:]) 

1001 raise _ellipticError(n, sn, cn, dn, cause=e) 

1002 

1003 

1004def _ellipticError(where, *args, **kwds_cause_txt): 

1005 '''(INTERNAL) Format an L{EllipticError}. 

1006 ''' 

1007 def _x_t_kwds(cause=None, txt=NN, **kwds): 

1008 return cause, txt, kwds 

1009 

1010 x, t, kwds = _x_t_kwds(**kwds_cause_txt) 

1011 

1012 n = _dunder_nameof(where, where) 

1013 n = _DOT_(Elliptic.__name__, n) 

1014 n = _SPACE_(_invokation_, n) 

1015 u = unstr(n, *args, **kwds) 

1016 return EllipticError(u, cause=x, txt=t) 

1017 

1018 

1019def _Horner(S, e1, E2, E3, E4, E5, *over): 

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

1021 ''' 

1022 E22 = E2**2 

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

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

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

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

1027 # converted to Horner-like form ... 

1028 e = e1 * 4084080 

1029 S *= e 

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

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

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

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

1034 S += 4084080 

1035 if over: 

1036 e *= over[0] 

1037 return S.fdiv(e) # Fsum 

1038 

1039 

1040def _iterations(inst, i): 

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

1042 ''' 

1043 if inst and i > 0: 

1044 inst._iteration += i 

1045 

1046 

1047def _3over(a, b): 

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

1049 ''' 

1050 return _over(_3_0, a * b) 

1051 

1052 

1053def _rC(unused, x, y): 

1054 '''(INTERNAL) Defined only for C{y != 0} and C{x >= 0}. 

1055 ''' 

1056 d = x - y 

1057 if d < 0: # catch NaN 

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

1059 d = -d 

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

1061 elif d == 0: # XXX d < EPS0? or EPS02 or _EPSmin 

1062 d, r = y, _1_0 

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

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

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

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

1067 else: # PYCHOK no cover 

1068 raise _ellipticError(Elliptic.fRC, x, y) 

1069 return r / sqrt(d) # float 

1070 

1071 

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

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

1074 ''' 

1075 L = _List(x, y, z) 

1076 S = _Dsum() 

1077 for a, m, r, s in L.amrs4(inst, True, _TolRF): 

1078 if s: 

1079 S += _over(_3_0, (r + z) * s[2] * m) 

1080 z = L[2] # s[2] = sqrt(z) 

1081 x, y = L.rescale(-a * m, x, y) 

1082 xy = x * y 

1083 z = (x + y) / _3_0 

1084 z2 = z**2 

1085 return _Horner(S(_1_0), sqrt(a) * a * m, 

1086 xy - _6_0 * z2, 

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

1088 (xy - z2) * _3_0 * z2, 

1089 xy * z2 * z, *over) # Fsum 

1090 

1091 

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

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

1094 ''' 

1095 for a, b in _ab2(inst, x, y): # PYCHOK yield 

1096 pass 

1097 return _over(PI, a + b) # float 

1098 

1099 

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

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

1102 ''' 

1103 L = _List(x, y, z) 

1104 for a, m, _, _ in L.amrs4(inst, False, _TolRF): 

1105 pass 

1106 x, y = L.rescale(a * m, x, y) 

1107 z = neg(x + y) 

1108 xy = x * y 

1109 e2 = xy - z**2 

1110 e3 = xy * z 

1111 e4 = e2**2 

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

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

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

1115 # converted to Horner-like form ... 

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

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

1118 S += 240240 

1119 return S.fdiv(sqrt(a) * 240240) # Fsum 

1120 

1121 

1122def _rG2(inst, x, y, PI_=PI_4): # 2-args 

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

1124 ''' 

1125 m = -1 # neg! 

1126 S = None 

1127 for a, b in _ab2(inst, x, y): # PYCHOK yield 

1128 if S is None: # initial 

1129 S = _Dsum() 

1130 S += (a + b)**2 * _0_5 

1131 else: 

1132 S += (a - b)**2 * m 

1133 m *= 2 

1134 return S(PI_).fover(a + b) 

1135 

1136 

1137def _rG3(inst, x, y, z): # 3-arg version 

1138 '''(INTERNAL) C{x}, C{y} and C{z} all non-zero, see C{.fRG}. 

1139 ''' 

1140 R = _RF3(inst, x, y, z) * z 

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

1142 if rd: # Carlson, eq 1.7 

1143 R += _RD(inst, x, y, z, _3_0 / rd) 

1144 R += sqrt(x * y / z) 

1145 return R.fover(_2_0) 

1146 

1147 

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

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

1150 ''' 

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

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

1153 

1154 L = _List(x, y, z, p) 

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

1156 S = _Dsum() 

1157 for a, m, _, s in L.amrs4(inst, True, _TolRD): 

1158 if s: 

1159 d = _xyzp(*s) 

1160 if d: 

1161 if n: 

1162 rc = _rC(inst, _1_0, n / d**2 + _1_0) 

1163 n *= _1_64th # /= chokes PyChecker 

1164 else: 

1165 rc = _1_0 # == _rC(None, _1_0, _1_0) 

1166 S += rc / (d * m) 

1167 else: # PYCHOK no cover 

1168 return NAN 

1169 x, y, z = L.rescale(a * m, x, y, z) 

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

1171 p2 = p**2 

1172 p3 = p2 * p 

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

1174 E2p = E2 * p 

1175 xyz = x * y * z 

1176 return _Horner(S(_6_0), sqrt(a) * a * m, E2, 

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

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

1179 xyz * p2, *over) # Fsum 

1180 

1181# **) MIT License 

1182# 

1183# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

1184# 

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

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

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

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

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

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

1191# 

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

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

1194# 

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

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

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

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

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

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

1201# OTHER DEALINGS IN THE SOFTWARE.