1

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

3

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

5

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

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

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

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

10

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

12

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

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

15

16 - member functions with 1 or 3 arguments are renamed to Python methods

17 starting with an C{"f"}, example C{E(psi)} to C{fE(psi)} and C{E(sn,

18 cn, dn)} to C{fE(sn, cn, dn)},

19

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

21 letter or an underscore if private.

22

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

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

25 file C{Header}.

26

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

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

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

30

31 B{Elliptic integrals and functions.}

32

33 This provides the elliptic functions and integrals needed for

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

35 categories of function are provided:

36

37 - functions to compute U{symmetric elliptic integrals

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

39

40 - methods to compute U{Legrendre's elliptic integrals

41 <https://DLMF.NIST.gov/19.2.ii>} and U{Jacobi elliptic

42 functions<https://DLMF.NIST.gov/22.2>}.

43

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

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

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

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

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

49

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

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

52

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

54

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

56 period I{C{pi}}.

57

58 The computation of the elliptic integrals uses the algorithms given

59 in U{B. C. Carlson, Computation of real or complex elliptic integrals

60<https://DOI.org/10.1007/BF02198293>} (also available U{here

61<https://ArXiv.org/pdf/math/9409227.pdf>}), Numerical Algorithms 10,

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

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

64

65 The computation of the Jacobi elliptic functions uses the algorithm

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

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

68 Numerische Mathematik 7, 78--90 (1965) or optionally the C{Jacobi

69 amplitude} in method L{Elliptic.sncndn<pygeodesy.Elliptic.sncndn>}.

70

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

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

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

74'''

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

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

77

78 from pygeodesy.basics import copysign0, map2, neg, neg_

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

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

81 _4_0, _6_0, _8_0, _64_0, _180_0, _360_0, \

82 _EPStol as _TolJAC, _over, \

83 _EPSjam as _TolJAM # PYCHOK used!

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

85 from pygeodesy.fmath import favg, hypot1, zqrt, _ValueError

86 from pygeodesy.fsums import Fsum, _sum

87 from pygeodesy.internals import _DUNDER_nameof

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

89 _invokation_, _negative_, _SPACE_

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

91# from pygeodesy.lazily import _ALL_LAZY # from .named

92 from pygeodesy.named import _Named, _NamedTuple, _ALL_LAZY, Fmt, unstr

93 from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all

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

95 from pygeodesy.units import Scalar, Scalar_

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

97

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

99 radians, sin, sinh, sqrt, tan, tanh

100

101__all__ = _ALL_LAZY.elliptic

102__version__ = '24.10.14'

103

104_TolRD = zqrt(EPS * 0.002)

105_TolRF = zqrt(EPS * 0.030)

106_TolRG0 = _TolJAC * 2.7

107_TRIPS = 25 # Max depth, 6-7 might be sufficient

108

109

110 class _Cs(object):

111 '''(INTERAL) Complete integrals cache.

112 '''

113 def __init__(self, **kwds):

114 self.__dict__ = kwds

115

116

117 class _Dsum(list):

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

119 '''

120 def __call__(self, s):

121 try: # Fsum *= s

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

123 except ValueError: # Fsum(NAN) exception

124 return _sum(self) * s

125

126 def __iadd__(self, x):

127 list.append(self, x)

128 return self

129

130

131 class _Dadd(_Dsum):

132 '''(INTERNAL) Deferred C{Fsum} for C{_List.amrs4}.

133 '''

134 def __init__(self, mul):

135 self._mul = mul

136

137 def __add__(self, x):

138 self += x

139 r = self(self._mul)

140 _ = self.pop()

141 return r # Fsum or float

142

143

144 class Elliptic(_Named):

145 '''Elliptic integrals and functions.

146

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

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

149 '''

150# _alpha2 = 0

151# _alphap2 = 0

152# _eps = EPS

153# _k2 = 0

154# _kp2 = 0

155

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

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

158

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

160

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

162

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

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

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

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

167 '''

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

169 if name:

170 self.name = name

171

172 @Property_RO

173 def alpha2(self):

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

175 '''

176 return self._alpha2

177

178 @Property_RO

179 def alphap2(self):

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

181 '''

182 return self._alphap2

183

184 @Property_RO

185 def cD(self):

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

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

188 '''

189 return self._cDEKEeps.cD

190

191 @Property_RO

192 def _cDEKEeps(self):

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

194 '''

195 k2, kp2 = self.k2, self.kp2

196 if k2:

197 if kp2:

198 try:

199 self._iteration = 0

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

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

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

203 cD = float(D)

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

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

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

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

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

209 cK = _rF2(self, kp2, _1_0)

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

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

212

213 except Exception as e:

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

215 else:

216 cD = cK = cKE = INF

217 cE = _1_0

218 eps = k2

219 else:

220 cD = PI_4

221 cE = cK = PI_2

222 cKE = _0_0 # k2 * cD

223 eps = EPS

224

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

226

227 @Property_RO

228 def cE(self):

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

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

231 '''

232 return self._cDEKEeps.cE

233

234 @Property_RO

235 def cG(self):

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

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

238 '''

239 return self._cGHPi.cG

240

241 @Property_RO

242 def _cGHPi(self):

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

244 '''

245 alpha2, alphap2, kp2 = self.alpha2, self.alphap2, self.kp2

246 try:

247 self._iteration = 0

248 if alpha2:

249 if alphap2:

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

251 cK = self.cK

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

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

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

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

256 else:

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

258 cPi = INF # XXX or NAN?

259 else:

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

261 else:

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

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

264 # So write (for alpha2 = 0)

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

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

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

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

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

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

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

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

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

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

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

276

277 except Exception as e:

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

279 alphap2=alphap2, cause=e)

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

281

282 @Property_RO

283 def cH(self):

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

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

286 '''

287 return self._cGHPi.cH

288

289 @Property_RO

290 def cK(self):

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

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

293 '''

294 return self._cDEKEeps.cK

295

296 @Property_RO

297 def cKE(self):

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

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

300 '''

301 return self._cDEKEeps.cKE

302

303 @Property_RO

304 def cPi(self):

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

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

307 '''

308 return self._cGHPi.cPi

309

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

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

312

313 @arg sn: sin(φ).

314 @arg cn: cos(φ).

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

316

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

318

319 @raise EllipticError: Invalid invokation or no convergence.

320 '''

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

322

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

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

325

326 @arg sn: sin(φ).

327 @arg cn: cos(φ).

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

329

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

331

332 @raise EllipticError: Invalid invokation or no convergence.

333 '''

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

335

336 def deltaEinv(self, stau, ctau):

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

338

339 @arg stau: sin(τ)

340 @arg ctau: cos(τ)

341

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

343

344 @raise EllipticError: No convergence.

345 '''

346 try:

347 if _signBit(ctau): # pi periodic

348 stau, ctau = neg_(stau, ctau)

349 t = atan2(stau, ctau)

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

351

352 except Exception as e:

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

354

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

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

357

358 @arg sn: sin(φ).

359 @arg cn: cos(φ).

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

361

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

363

364 @raise EllipticError: Invalid invokation or no convergence.

365 '''

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

367

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

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

370

371 @arg sn: sin(φ).

372 @arg cn: cos(φ).

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

374

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

376

377 @raise EllipticError: Invalid invokation or no convergence.

378 '''

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

380

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

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

383

384 @arg sn: sin(φ).

385 @arg cn: cos(φ).

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

387

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

389

390 @raise EllipticError: Invalid invokation or no convergence.

391 '''

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

393

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

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

396

397 @arg sn: sin(φ).

398 @arg cn: cos(φ).

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

400

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

402 (C{float}).

403

404 @raise EllipticError: Invalid invokation or no convergence.

405 '''

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

407

408 def _Einv(self, x):

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

410 '''

411 E2 = self.cE * _2_0

412 n = floor(x / E2 + _0_5)

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

414 # linear approximation

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

416 Phi = Fsum(phi)

417 # first order correction

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

419 self._iteration = 0

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

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

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

423 _Phi2 = Phi.fsum2f_ # aggregate

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

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

426 if dn:

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

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

429 else: # PYCHOK no cover

430 d = _0_0 # XXX or continue?

431 if fabs(d) < _TolJAC: # 3-4 trips

432 _iterations(self, i)

433 break

434 else: # PYCHOK no cover

435 raise _convergenceError(d, _TolJAC)

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

437

438 @Property_RO

439 def eps(self):

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

441 '''

442 return self._cDEKEeps.eps

443

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

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

446 Jacobi elliptic functions.

447

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

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

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

451

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

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

454

455 @raise EllipticError: Invalid invokation or no convergence.

456 '''

457 def _fD(sn, cn, dn):

458 r = fabs(sn)**3

459 if r:

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

461 return r

462

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

464 self.deltaD, _fD)

465

466 def fDelta(self, sn, cn):

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

468

469 @arg sn: sin(φ).

470 @arg cn: cos(φ).

471

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

473 '''

474 try:

475 k2 = self.k2

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

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

478 return sqrt(s) if s else _0_0

479

480 except Exception as e:

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

482

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

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

485 Jacobi elliptic functions.

486

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

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

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

490

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

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

493

494 @raise EllipticError: Invalid invokation or no convergence.

495 '''

496 def _fE(sn, cn, dn):

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

498 '''

499 if sn:

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

501 kp2, k2 = self.kp2, self.k2

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

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

504 if k2:

505 Ei -= _RD(self, cn2, dn2, _1_0, _3over(k2, sn**2))

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

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

508 if kp2:

509 Ei += (_RD( self, cn2, _1_0, dn2, _3over(k2, sn**2)) +

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

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

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

513 Ei -= _RD(self, dn2, _1_0, cn2, _3over(kp2, sn**2))

514 Ei *= fabs(sn)

515 ei = float(Ei)

516 else: # PYCHOK no cover

517 ei = _0_0

518 return ei

519

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

521 self.deltaE, _fE,

522 k2_0=self.k2==0)

523

524 def fEd(self, deg):

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

526 the argument given in C{degrees}.

527

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

529

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

531

532 @raise EllipticError: No convergence.

533 '''

534 if _K_2_0:

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

536 elif fabs(deg) < _180_0:

537 e = _0_0

538 else:

539 e = ceil(deg / _360_0 - _0_5)

540 deg -= e * _360_0

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

542

543 def fEinv(self, x):

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

545

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

547

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

549 (C{float}).

550

551 @raise EllipticError: No convergence.

552 '''

553 try:

554 return self._Einv(x)

555 except Exception as e:

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

557

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

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

560 Jacobi elliptic functions.

561

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

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

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

565

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

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

568

569 @raise EllipticError: Invalid invokation or no convergence.

570 '''

571 def _fF(sn, cn, dn):

572 r = fabs(sn)

573 if r:

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

575 return r

576

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

578 self.deltaF, _fF,

579 k2_0=self.k2==0, kp2_0=self.kp2==0)

580

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

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

583 Jacobi elliptic functions.

584

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

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

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

588

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

590

591 @raise EllipticError: Invalid invokation or no convergence.

592

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

594 geodesic in terms of this combination of elliptic

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

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

597 riIOAAAAQAAJ&pg=PA181>}.

598

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

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

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

602 '''

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

604 self.cG, self.deltaG)

605

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

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

608 Jacobi elliptic functions.

609

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

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

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

613

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

615

616 @raise EllipticError: Invalid invokation or no convergence.

617

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

619 on the geodesic in terms of this combination of

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

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

622

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

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

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

626 '''

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

628 self.cH, self.deltaH)

629

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

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

632 Jacobi elliptic functions.

633

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

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

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

637

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

639

640 @raise EllipticError: Invalid invokation or no convergence.

641 '''

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

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

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

645 self.cPi, self.deltaPi)

646

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

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

649 '''

650 def _fX(sn, cn, dn):

651 if sn:

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

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

654 if aX:

655 sn2 = sn**2

656 p = sn2 * self.alphap2 + cn2

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

658 R *= fabs(sn)

659 r = float(R)

660 else: # PYCHOK no cover

661 r = _0_0

662 return r

663

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

665

666 def _fXf(self, phi_or_sn, cn, dn, cX, deltaX, fX, k2_0=False, kp2_0=False):

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

668 '''

669 self._iteration = 0 # aggregate

670 phi = sn = phi_or_sn

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

672 if k2_0: # C++ version 2.4

673 return phi

674 elif kp2_0:

675 return asinh(tan(phi))

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

677 if fabs(phi) >= PI:

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

679 # fall through

680 elif cn is None or dn is None:

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

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

683

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

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

686 else:

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

688 return copysign0(xi, sn)

689

690 def _jam(self, x):

691 '''Jacobi amplitude function.

692

693 @return: C{phi} (C{float}).

694 '''

695 # implements DLMF Sec 22.20(ii), see also U{Sala

696 # (1989)<https://doi.org/10.1137/0520100>}, Sec 5

697 if self.k2:

698 if self.kp2:

699 r, ac = self._jamac2

700 x *= r # phi

701 for a, c in ac:

702 p = x

703 x = favg(asin(c * sin(x) / a), x)

704 if self.k2 < 0: # Sala Eq. 5.8

705 x = p - x

706 else: # PYCHOK no cover

707 x = atan(sinh(x)) # gd(x)

708 return x

709

710 @Property_RO

711 def _jamac2(self):

712 '''(INTERNAL) Get Jacobi amplitude 2-tuple C{(r, ac)}.

713 '''

714 a = r = _1_0

715 b, c = self.kp2, self.k2

716 # assert b and c

717 if c < 0: # Sala Eq. 5.8

718 r = sqrt(b)

719 b = _1_0 / b

720# c *= b # unused

721 ac = [] # [(a, sqrt(c))] unused

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

723 b = sqrt(a * b)

724 c = favg(a, -b)

725 a = favg(a, b) # == PI_2 / K

726 ac.append((a, c))

727 if c <= (a * _TolJAM): # 7 trips, quadratic

728 _iterations(self, i)

729 break

730 else: # PYCHOK no cover

731 raise _convergenceError(c / a, _TolJAM)

732 r *= a * float(1 << i) # 2**i == 2**len(ac)

733 return r, tuple(reversed(ac))

734

735 @Property_RO

736 def k2(self):

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

738 '''

739 return self._k2

740

741 @Property_RO

742 def kp2(self):

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

744 '''

745 return self._kp2

746

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

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

749

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

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

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

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

754

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

756 or B{C{alphap2}}.

757

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

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

760 that these conditions are met to enable accuracy to be

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

762 '''

763 if self.__dict__:

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

765

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

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

768

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

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

771 Error=EllipticError)

772

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

774 # K E D

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

776 # k = 1: inf 1 inf

777 # Pi G H

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

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

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

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

782 #

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

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

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

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

787 # Pi(alpha2, 1) = inf

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

789

790 def sncndn(self, x, jam=False):

791 '''The Jacobi amplitude and elliptic function.

792

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

794 @kwarg jam: If C{True}, use the Jacobi amplitude otherwise the

795 Bulirsch' function (C{bool}).

796

797 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with C{*n(B{x}, k)}.

798

799 @raise EllipticError: No convergence.

800 '''

801 self._iteration = 0 # reset

802 try:

803 if self.kp2:

804 if jam: # Jacobi amplitude, C++ version 2.4

805 sn, cn, dn = self._sncndn3(self._jam(x))

806

807 else: # Bulirsch's sncndn routine, p 89 of

808 # Numerische Mathematik 7, 78-90 (1965).

809 # Implements DLMF Eqs 22.17.2 - 22.17.4,

810 # but only good for .k2 > 1 or .kp2 < 0

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

812 dn = _1_0

813 sn, cn = _sincos2(x * cd)

814 if sn:

815 a = cn / sn

816 c *= a

817 for m, n in mn:

818 a *= c

819 c *= dn

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

821 a = c / m

822 a = _1_0 / hypot1(c)

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

824 cn = sn * c

825 if d: # and _signBit(self.kp2): # implied

826 cn, dn = dn, cn

827 sn = sn / d # /= chokes PyChecker

828 else:

829 sn = tanh(x) # accurate for large abs(x)

830 cn = dn = _1_0 / cosh(x)

831

832 except Exception as e:

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

834

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

836

837 def _sncndn3(self, phi):

838 '''(INTERNAL) Helper for C{.fEinv}, C{._fXf} and C{.sncndn}.

839 '''

840 sn, cn = _sincos2(phi)

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

842

843 @Property_RO

844 def _sncndn4(self):

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

846 '''

847 d, b = 0, self.kp2

848 if _signBit(b):

849 d = _1_0 - b

850 b = neg(b / d)

851 d = sqrt(d)

852 ab, a = [], _1_0

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

854 b = sqrt(b)

855 ab.append((a, b))

856 c = favg(a, b)

857 r = fabs(a - b)

858 if r <= (a * _TolJAC): # 6 trips, quadratic

859 _iterations(self, i)

860 break

861 t = a

862 b *= a

863 a = c

864 else: # PYCHOK no cover

865 raise _convergenceError(r / t, _TolJAC)

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

867 return c, d, cd, tuple(reversed(ab))

868

869 @staticmethod

870 def fRC(x, y):

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

872

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

874

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

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

877 '''

878 return _rC(None, x, y)

879

880 @staticmethod

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

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

883

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

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

886

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

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

889 '''

890 try:

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

892 except Exception as e:

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

894

895 @staticmethod

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

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

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

899

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

901

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

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

904 '''

905 try:

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

907 except Exception as e:

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

909

910 @staticmethod

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

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

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

914

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

916

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

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

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

920 EllipticFunction_8cpp_source.html#l00096>} version 2.3.

921 '''

922 try:

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

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

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

926 except Exception as e:

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

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

929

930 @staticmethod

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

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

933

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

935

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

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

938 '''

939 try:

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

941 except Exception as e:

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

943

944 @staticmethod

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

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

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

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

949 if m:

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

951 except Exception as e:

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

953 return float(R)

954

955_allPropertiesOf_n(16, Elliptic) # # PYCHOK assert, see Elliptic.reset

956

957

958 class EllipticError(_ValueError):

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

960 '''

961 pass

962

963

964 class Elliptic3Tuple(_NamedTuple):

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

966 '''

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

968 _Units_ = ( Scalar, Scalar, Scalar)

969

970

971 class _List(list):

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

973 '''

974 _a0 = None

975

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

977 list.__init__(self, xyzp)

978

979 def a0(self, n):

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

981 '''

982 A = _Dsum(self)

983 while len(A) < n:

984 A += A[-1]

985 self._a0 = a = float(A(_1_0 / n))

986 return a

987

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

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

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

991 '''

992 L = self

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

994 m = 1

995 t = max(fabs(a - _) for _ in L) / Tol # thresh

996 for i in range(_TRIPS):

997 d = fabs(a * m)

998 if d > t: # 3-6 trips

999 _iterations(inst, i)

1000 break

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

1002 # Deferred fdot(s[:3], s[1], s[2], s[0]) + ...

1003 Q = _Dadd(_0_25)

1004 Q += s[0] * s[1]

1005 Q += s[1] * s[2]

1006 Q += s[2] * s[0]

1007 L[:] = (float(Q + _) for _ in L)

1008 a = float(Q + a)

1009 if y: # yield only if used

1010 r = _sum(Q) # fdot

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

1012 m *= 4

1013 else: # PYCHOK no cover

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

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

1016

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

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

1019 '''

1020 # assert am

1021 a0 = self._a0

1022 _am = _1_0 / am

1023 for x in xs:

1024 yield (a0 - x) * _am

1025

1026

1027 def _ab2(inst, x, y):

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

1029 '''

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

1031 if b > a:

1032 b, a = a, b

1033 for i in range(_TRIPS):

1034 yield a, b # xi, yi

1035 d = fabs(a - b)

1036 if d <= (a * _TolRG0): # 3-4 trips

1037 _iterations(inst, i)

1038 break

1039 t = a

1040 a = favg(t, b)

1041 b = sqrt(t * b)

1042 else: # PYCHOK no cover

1043 raise _convergenceError(d / t, _TolRG0)

1044

1045

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

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

1048 '''

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

1050 return ValueError(t) # txt only

1051

1052

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

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

1055 '''

1056 try:

1057 if cn is None or dn is None:

1058 raise ValueError(_invalid_)

1059

1060 if _signBit(cn):

1061 sn, cn = neg_(sn, cn)

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

1063 return r - atan2(sn, cn)

1064

1065 except Exception as e:

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

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

1068

1069

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

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

1072 '''

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

1074 return cause, txt, kwds

1075

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

1077

1078 n = _DUNDER_nameof(where, where)

1079 n = _DOT_(Elliptic.__name__, n)

1080 n = _SPACE_(_invokation_, n)

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

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

1083

1084

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

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

1087 '''

1088 E22 = E2**2

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

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

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

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

1093 # converted to Horner-like form ...

1094 e = e1 * 4084080

1095 S *= e

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

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

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

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

1100 S += 4084080

1101 if over:

1102 e *= over[0]

1103 return S.fdiv(e) # Fsum

1104

1105

1106 def _iterations(inst, i):

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

1108 '''

1109 if inst and i > 0:

1110 inst._iteration += i

1111

1112

1113 def _3over(a, b):

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

1115 '''

1116 return _over(_3_0, a * b)

1117

1118

1119 def _rC(unused, x, y):

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

1121 '''

1122 d = x - y

1123 if d < 0: # catch NaN

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

1125 d = -d

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

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

1128 d, r = y, _1_0

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

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

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

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

1133 else: # PYCHOK no cover

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

1135 return r / sqrt(d) # float

1136

1137

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

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

1140 '''

1141 L = _List(x, y, z)

1142 S = _Dsum()

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

1144 if s:

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

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

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

1148 xy = x * y

1149 z = (x + y) / _3_0

1150 z2 = z**2

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

1152 (xy - z2 * _6_0),

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

1154 (xy - z2) * z2 * _3_0,

1155 (xy * z2 * z), *over) # Fsum

1156

1157

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

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

1160 '''

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

1162 pass

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

1164

1165

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

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

1168 '''

1169 L = _List(x, y, z)

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

1171 pass

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

1173 z = neg(x + y)

1174 xy = x * y

1175 e2 = xy - z**2

1176 e3 = xy * z

1177 e4 = e2**2

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

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

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

1181 # converted to Horner-like form ...

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

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

1184 S += 240240

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

1186

1187

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

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

1190 '''

1191 m = -1 # neg!

1192 S = None

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

1194 if S is None: # initial

1195 S = _Dsum()

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

1197 else:

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

1199 m *= 2

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

1201

1202

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

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

1205 '''

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

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

1208 if rd: # Carlson, eq 1.7

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

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

1211 return R.fover(_2_0)

1212

1213

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

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

1216 '''

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

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

1219

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

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

1222 S = _Dsum()

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

1224 if s:

1225 d = _xyzp(*s)

1226 if d:

1227 if n:

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

1229 n = n / _64_0 # /= chokes PyChecker

1230 else:

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

1232 S += rc / (d * m)

1233 else: # PYCHOK no cover

1234 return NAN

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

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

1237 p2 = p**2

1238 p3 = p2 * p

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

1240 E2p = E2 * p

1241 xyz = x * y * z

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

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

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

1245 xyz * p2, *over) # Fsum

1246

1247# **) MIT License

1248#

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

1250#

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

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

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

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

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

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

1257#

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

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

1260#

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

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

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

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

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

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

1267# OTHER DEALINGS IN THE SOFTWARE.