Coverage for pygeodesy/elliptic.py: 97%

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

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

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

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

8including symmetric integrals L{Elliptic.fRC}, L{Elliptic.fRD}, L{Elliptic.fRF},

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

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

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

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

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

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

19 - other Python method names conventionally start with a lower-case letter or an

20 underscore if private.

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

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

24file C{Header}.

26Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2024) and licensed

27under the MIT/X11 License. For more information, see the U{GeographicLib

28<https://GeographicLib.SourceForge.io>} documentation.

30B{Elliptic integrals and functions.}

32This provides the elliptic functions and integrals needed for classes C{Ellipsoid},

33C{GeodesicExact}, and C{TransverseMercatorExact}. Two categories of functions are

34provided:

36 - functions to compute U{symmetric elliptic integrals<https://DLMF.NIST.gov/19.16.i>}

38 - methods to compute U{Legrendre's elliptic integrals<https://DLMF.NIST.gov/19.2.ii>}

39 and U{Jacobi elliptic functions<https://DLMF.NIST.gov/22.2>}.

41In the latter case, an object is constructed giving the modulus C{k} (and optionally

42the parameter C{alpha}). The modulus (and parameter) are always passed as squares

43which allows C{k} to be pure imaginary. (Confusingly, Abramowitz and Stegun call

44C{m = k**2} the "parameter" and C{n = alpha**2} the "characteristic".)

46In geodesic applications, it is convenient to separate the incomplete integrals into

47secular and periodic components, e.g.

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

51where I{C{delta E(phi, k)}} is an odd periodic function with period I{C{pi}}.

53The computation of the elliptic integrals uses the algorithms given in U{B. C. Carlson,

54Computation of real or complex elliptic integrals <https://DOI.org/10.1007/BF02198293>}

55(also available U{here<https://ArXiv.org/pdf/math/9409227.pdf>}), Numerical Algorithms

5610, 13--26 (1995) with the additional optimizations given U{here<https://DLMF.NIST.gov/19.36.i>}.

58The computation of the Jacobi elliptic functions uses the algorithm given in U{R. Bulirsch, Numerical

59Calculation of Elliptic Integrals and Elliptic Functions<https://DOI.org/10.1007/BF01397975>},

60Numerische Mathematik 7, 78--90 (1965) or optionally the C{Jacobi amplitude} in method

61L{Elliptic.sncndn<pygeodesy.Elliptic.sncndn>}.

63The notation follows U{NIST Digital Library of Mathematical Functions <https://DLMF.NIST.gov>}

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

65'''

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

67from __future__ import division as _; del _ # noqa: E702 ;

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

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

71 _0_0, _0_25, _0_5, _1_0, _2_0, _3_0, _4_0, _6_0, _8_0, \

72 _64_0, _180_0, _360_0, _flipsign, _over, _1_over

73from pygeodesy.constants import _EPSjam as _TolJAM # PYCHOK used!

74# from pygeodesy.ellipsoids import Ellipsoid # _MODS

75from pygeodesy.errors import _ConvergenceError, _ValueError

76from pygeodesy.fmath import favg, Fdot, Fdot_, hypot1, zqrt, _operator

77from pygeodesy.fsums import Fsum, _fsum

78from pygeodesy.internals import _Enum, typename

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

80 _negative_, _SPACE_

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

82# from pygeodesy.lazily import _ALL_LAZY # from .utily

83from pygeodesy.named import _Named, _NamedTuple, unstr

84from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all

85# from pygeodesy.streprs import unstr # from .named

86# from pygeodesy.triaxials import Triaxial_ # _MODS

87from pygeodesy.units import Scalar, Scalar_

88from pygeodesy.utily import atan2, _ALL_LAZY # sincos2 as _sincos2

90from math import asin, asinh, atan, ceil, cosh, fabs, floor, radians, \

91 sin, sinh, sqrt, tan, tanh # tan as _tan

92# import operator as _operator # from .fmath

94__all__ = _ALL_LAZY.elliptic

95__version__ = '26.02.15'

97_TolRD = zqrt(EPS * 0.002)

98_TolRF = zqrt(EPS * 0.030)

99_TolRG0 = _TolJAC * 2.7

100_MAXIT = 32 # ._B4 max depth, 6..18 sufficient

101

102

103class _Cs(_Enum):

104 '''(INTERAL) Complete Integrals cache.

105 '''

106 pass

107

108

109class Elliptic(_Named):

110 '''Elliptic integrals and functions.

111

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

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

114 '''

115# _alpha2 = 0

116# _alphap2 = 0

117# _eps = EPS

118# _k2 = 0

119# _kp2 = 0

120

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

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

123

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

125

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

127

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

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

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

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

132 '''

133 if name:

134 self.name = name

135 self._reset(k2, alpha2, kp2, alphap2)

136

137 @Property_RO

138 def alpha2(self):

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

140 '''

141 return self._alpha2

142

143 @Property_RO

144 def alphap2(self):

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

146 '''

147 return self._alphap2

148

149 def _B3(self, x):

150 '''(INTERNAL) Bulirsch' sncndn routine.

151 '''

152 # Numerische Mathematik 7, 1965, P 89

153 # implements DLMF Eqs 22.17.2 - 22.17.4

154 c, d, cd, mn = self._B4

155 sn, cn = _sincos2(x * cd)

156 dn = _1_0

157 if sn:

158 a = cn / sn

159 c *= a

160 for m, n in mn:

161 a *= c

162 c *= dn

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

164 a = c / m

165 a = _1_over(hypot1(c))

166 sn = _flipsign(a, sn)

167 cn = sn * c

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

169 cn, dn = dn, cn

170 sn = sn / d # /= chokes PyChecker

171 return sn, cn, dn

172

173 @Property_RO

174 def _B4(self):

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

176 '''

177 a = P = _1_0

178 b, d = self.kp2, 0 # kp2 >= 0 always here

179 if _signBit(b): # PYCHOK no cover

180 d = a - b

181 b = neg(b / d)

182 d = sqrt(d)

183 mn = []

184 _mn = mn.append

185 for i in range(1, _MAXIT): # GEOGRAPHICLIB_PANIC

186 b = sqrt(P * b)

187 _mn((a, b))

188 c = favg(a, b)

189 r = fabs(a - b)

190 if r <= (a * _TolJAC): # 4..6 trips, quadratic

191 self._iteration += i

192 break

193 P, a = a, c

194 else: # PYCHOK no cover

195 raise _ConvergenceError(_MAXIT, _over(r, P), _TolJAC)

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

197 return c, d, cd, tuple(reversed(mn))

198

199 @Property_RO

200 def cD(self):

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

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

203 '''

204 return self._cDEKEeps.cD

205

206 @Property_RO

207 def _cDEKEeps(self):

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

209 '''

210 k2, kp2 = self.k2, self.kp2

211 if k2:

212 if kp2:

213 try:

214 # self._iteration = 0

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

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

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

218 cD = float(D)

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

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

221 cE = self._cE(kp2)

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

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

224 cK = self._cK(kp2)

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

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

227

228 except Exception as X:

229 raise _ellipticError(self, k2=k2, kp2=kp2, cause=X)

230 else:

231 cD = cK = cKE = INF

232 cE = _1_0

233 eps = k2

234 else:

235 cD = PI_4

236 cE = cK = PI_2

237 cKE = _0_0 # k2 * cD

238 eps = EPS

239

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

241

242 @Property_RO

243 def cE(self):

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

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

246 '''

247 return self._cDEKEeps.cE

248

249 def _cE(self, kp2): # compl integr 2nd kind

250 return _rG2(kp2, _1_0, self, PI_=PI_2)

251

252 @Property_RO

253 def cG(self):

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

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

256 '''

257 return self._cGHPi.cG

258

259 @Property_RO

260 def _cGHPi(self):

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

262 '''

263 alpha2, alphap2, kp2 = self.alpha2, self.alphap2, self.kp2

264 try:

265 # self._iteration = 0

266 if alpha2:

267 if alphap2:

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

269 cK = self.cK

270 Rj = _RJfma(_0_0, kp2, _1_0, alphap2, _3_0, self)

271 cG = Rj.ma(alpha2 - self.k2, cK) # G(alpha2, k)

272 cH = -Rj.ma(alphap2, -cK) # H(alpha2, k)

273 cPi = Rj.ma(alpha2, cK) # Pi(alpha2, k)

274 else:

275 cG = cH = _rC(_1_0, alphap2)

276 cPi = INF # XXX or NAN?

277 else:

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

279 else:

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

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

282 # So write (for alpha2 = 0)

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

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

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

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

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

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

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

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

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

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

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

294

295 except Exception as X:

296 raise _ellipticError(self, kp2=kp2, alpha2 =alpha2,

297 alphap2=alphap2, cause=X)

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

299

300 @Property_RO

301 def cH(self):

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

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

304 '''

305 return self._cGHPi.cH

306

307 @Property_RO

308 def cK(self):

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

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

311 '''

312 return self._cDEKEeps.cK

313

314 def _cK(self, kp2): # compl integr 1st kind

315 return _rF2(kp2, _1_0, self)

316

317 @Property_RO

318 def cKE(self):

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

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

321 '''

322 return self._cDEKEeps.cKE

323

324 @Property_RO

325 def cPi(self):

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

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

328 '''

329 return self._cGHPi.cPi

330

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

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

333

334 @arg sn: sin(φ).

335 @arg cn: cos(φ).

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

337

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

339

340 @raise EllipticError: Invalid invokation or no convergence.

341 '''

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

343

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

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

346

347 @arg sn: sin(φ).

348 @arg cn: cos(φ).

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

350

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

352

353 @raise EllipticError: Invalid invokation or no convergence.

354 '''

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

356

357 def deltaEinv(self, stau, ctau):

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

359

360 @arg stau: sin(τ)

361 @arg ctau: cos(τ)

362

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

364

365 @raise EllipticError: No convergence.

366 '''

367 try:

368 if _signBit(ctau): # pi periodic

369 stau, ctau = neg_(stau, ctau)

370 t = atan2(stau, ctau)

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

372

373 except Exception as X:

374 raise _ellipticError(self.deltaEinv, stau, ctau, cause=X)

375

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

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

378

379 @arg sn: sin(φ).

380 @arg cn: cos(φ).

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

382

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

384

385 @raise EllipticError: Invalid invokation or no convergence.

386 '''

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

388

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

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

391

392 @arg sn: sin(φ).

393 @arg cn: cos(φ).

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

395

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

397

398 @raise EllipticError: Invalid invokation or no convergence.

399 '''

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

401

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

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

404

405 @arg sn: sin(φ).

406 @arg cn: cos(φ).

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

408

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

410

411 @raise EllipticError: Invalid invokation or no convergence.

412 '''

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

414

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

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

417

418 @arg sn: sin(φ).

419 @arg cn: cos(φ).

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

421

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

423 (C{float}).

424

425 @raise EllipticError: Invalid invokation or no convergence.

426 '''

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

428

429 def _Einv(self, x):

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

431 '''

432 E2 = self.cE * _2_0

433 n = floor(x / E2 + _0_5)

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

435 # linear approximation

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

437 Phi = Fsum(phi)

438 # first order correction

439 phi = Phi.fsum_(-sin(phi * _2_0) * self.eps * _0_5)

440 # self._iteration = 0

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

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

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

444 _Phi2 = Phi.fsum2f_ # aggregate

445 for i in range(1, _MAXIT): # GEOGRAPHICLIB_PANIC

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

447 if dn:

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

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

450 else: # PYCHOK no cover

451 d = _0_0 # XXX continue?

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

453 self._iteration += i

454 break

455 else: # PYCHOK no cover

456 raise _ConvergenceError(_MAXIT, d, _TolJAC)

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

458

459 @Property_RO

460 def eps(self):

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

462 '''

463 return self._cDEKEeps.eps

464

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

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

467 Jacobi elliptic functions.

468

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

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

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

472

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

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

475

476 @raise EllipticError: Invalid invokation or no convergence.

477 '''

478 def _fD(sn, cn, dn):

479 r = fabs(sn)**3

480 if r:

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

482 return r

483

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

485 self.deltaD, _fD)

486

487 def fDelta(self, sn, cn):

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

489

490 @arg sn: sin(φ).

491 @arg cn: cos(φ).

492

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

494 '''

495 try:

496 k2 = self.k2

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

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

499 return sqrt(s) if s else _0_0

500

501 except Exception as X:

502 raise _ellipticError(self.fDelta, sn, cn, k2=k2, cause=X)

503

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

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

506 Jacobi elliptic functions.

507

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

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

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

511

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

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

514

515 @raise EllipticError: Invalid invokation or no convergence.

516 '''

517 def _fE(sn, cn, dn):

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

519 '''

520 if sn:

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

522 kp2, k2 = self.kp2, self.k2

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

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

525 if k2:

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

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

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

529 if kp2:

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

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

532 else: # PYCHOK no cover

533 Ei = _over(dn, fabs(cn)) # <https://DLMF.NIST.gov/19.25.E11>

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

535 Ei *= fabs(sn)

536 else: # PYCHOK no cover

537 Ei = _0_0

538 return float(Ei)

539

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

541 self.deltaE, _fE,

542 k2_0=self.k2==0)

543

544 def fEd(self, deg):

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

546 the argument given in C{degrees}.

547

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

549

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

551

552 @raise EllipticError: No convergence.

553 '''

554 if _K_2_0:

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

556 elif fabs(deg) < _180_0:

557 e = _0_0

558 else:

559 e = ceil(deg / _360_0 - _0_5)

560 deg -= e * _360_0

561 e *= self.cE * _4_0

562 return self.fE(radians(deg)) + e

563

564 def fEinv(self, x):

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

566

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

568

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

570 (C{float}).

571

572 @raise EllipticError: No convergence.

573 '''

574 try:

575 return self._Einv(x)

576 except Exception as X:

577 raise _ellipticError(self.fEinv, x, cause=X)

578

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

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

581 Jacobi elliptic functions.

582

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

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

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

586

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

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

589

590 @raise EllipticError: Invalid invokation or no convergence.

591 '''

592 def _fF(sn, cn, dn):

593 r = fabs(sn)

594 if r:

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

596 return r

597

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

599 self.deltaF, _fF,

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

601

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

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

604 Jacobi elliptic functions.

605

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

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

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

609

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

611

612 @raise EllipticError: Invalid invokation or no convergence.

613

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

615 geodesic in terms of this combination of elliptic

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

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

618 riIOAAAAQAAJ&pg=PA181>}.

619

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

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

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

623 '''

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

625 self.cG, self.deltaG)

626

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

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

629 Jacobi elliptic functions.

630

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

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

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

634

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

636

637 @raise EllipticError: Invalid invokation or no convergence.

638

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

640 on the geodesic in terms of this combination of

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

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

643

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

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

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

647 '''

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

649 self.cH, self.deltaH)

650

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

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

653 Jacobi elliptic functions.

654

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

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

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

658

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

660

661 @raise EllipticError: Invalid invokation or no convergence.

662 '''

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

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

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

666 self.cPi, self.deltaPi)

667

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

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

670 '''

671 def _fX(sn, cn, dn):

672 if sn:

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

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

675 if aX:

676 sn2 = sn**2

677 P = sn2 * self.alphap2 + cn2

678 R += _RJ(cn2, dn2, _1_0, P, _3over(aX, sn2), self)

679 R *= fabs(sn)

680 else: # PYCHOK no cover

681 R = _0_0

682 return float(R)

683

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

685

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

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

688 '''

689 # self._iteration = 0 # aggregate

690 phi = sn = phi_or_sn

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

692 if k2_0: # C++ version 2.4

693 return phi

694 elif kp2_0:

695 return asinh(tan(phi))

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

697 if fabs(phi) >= PI:

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

699 # fall through

700 elif cn is None or dn is None:

701 n = NN(_f_, typename(deltaX)[5:])

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

703

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

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

706 else:

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

708 return copysign0(xi, sn)

709

710 def _jam(self, x):

711 '''Jacobi amplitude function.

712

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

714 '''

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

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

717 if self.k2:

718 if self.kp2:

719 r, ac = self._jamac2

720 x *= r # phi

721 for a, c in ac:

722 P = x

723 a = asin(c * sin(x) / a)

724 x = favg(a, x)

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

726 x = P - x

727 else: # PYCHOK no cover

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

729 return x

730

731 @Property_RO

732 def _jamac2(self):

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

734 '''

735 a = r = _1_0

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

737 # assert b and c

738 if c < 0: # Sala Eq. 5.8

739 r = sqrt(b)

740 b = _1_over(b)

741# c *= b # unused

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

743 _ac = ac.append

744 for _ in range(_MAXIT): # GEOGRAPHICLIB_PANIC

745 b = sqrt(a * b)

746 c = favg(a, -b)

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

748 _ac((a, c))

749 if c <= (a * _TolJAM): # 7..18 trips, quadratic

750 break

751 else: # PYCHOK no cover

752 raise _ConvergenceError(_MAXIT, _over(c, a), _TolJAM)

753 i = len(ac)

754 r *= 2**i * a

755 self._iteration += i

756 return r, tuple(reversed(ac))

757

758 @Property_RO

759 def k2(self):

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

761 '''

762 return self._k2

763

764 @Property_RO

765 def kp2(self):

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

767 '''

768 return self._kp2

769

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

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

772

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

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

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

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

777

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

779 or B{C{alphap2}}.

780

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

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

783 that these conditions are met to enable accuracy to be

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

785 '''

786 _update_all(self, Base=Property_RO, needed=4)

787 self._reset(k2, alpha2, kp2, alphap2)

788

789 def _reset(self, k2, alpha2, kp2, alphap2):

790 '''(INITERNAL) Reset this elliptic.

791 '''

792 def _1p2(kp2, k2):

793 return (_1_0 - k2) if kp2 is None else kp2

794

795 def _S(**kwds):

796 return Scalar_(Error=EllipticError, **kwds)

797

798 self._k2 = _S(k2 = k2, low=None, high=_1_0)

799 self._kp2 = _S(kp2=_1p2(kp2, k2)) # low=_0_0

800

801 self._alpha2 = _S(alpha2 = alpha2, low=None, high=_1_0)

802 self._alphap2 = _S(alphap2=_1p2(alphap2, alpha2)) # low=_0_0

803

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

805 # K E D

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

807 # k = 1: inf 1 inf

808 # Pi G H

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

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

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

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

813 #

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

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

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

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

818 # Pi(alpha2, 1) = inf

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

820

821 self._iteration = 0

822

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

824 '''The Jacobi amplitude and elliptic function.

825

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

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

828 Bulirsch' function (C{bool}).

829

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

831

832 @raise EllipticError: No convergence.

833 '''

834 i = self._iteration

835 try:

836 if self.kp2:

837 if jam: # Jacobi amplitude, C++ v 2.4

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

839 else:

840 sn, cn, dn = self._B3(x)

841 else:

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

843 cn = dn = _1_over(cosh(x))

844

845 except Exception as X:

846 raise _ellipticError(self.sncndn, x, kp2=self.kp2, cause=X)

847

848 return Elliptic3Tuple(sn, cn, dn, iteration=self._iteration - i)

849

850 def _sncndn3(self, phi):

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

852 '''

853 sn, cn = _sincos2(phi)

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

855

856 @staticmethod

857 def fRC(x, y):

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

859

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

861

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

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

864 '''

865 return _rC(x, y)

866

867 @staticmethod

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

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

870

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

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

873

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

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

876 '''

877 try:

878 return float(_RD(x, y, z, *over))

879 except Exception as X:

880 raise _ellipticError(Elliptic.fRD, x, y, z, *over, cause=X)

881

882 @staticmethod

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

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

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

886

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

888

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

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

891 '''

892 try:

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

894 except Exception as X:

895 raise _ellipticError(Elliptic.fRF, x, y, z, cause=X)

896

897 @staticmethod

898 def _fRF3RD(x, z, m): # in .auxilats.AuxDLat.DE, -.AuxLat.Rectifying

899 y = _1_0 - m

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

901 R = _RF3(x, y, z)

902 if m:

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

904 except Exception as X:

905 raise _ellipticError(Elliptic._fRF3RD, x, y, z, m, cause=X)

906 return float(R)

907

908 @staticmethod

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

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

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

912

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

914

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

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

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

918 EllipticFunction_8cpp_source.html#l00096>} version 2.3.

919 '''

920 try:

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

922 _rG2(z, x) if y == 0 else (

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

924 except Exception as X:

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

926 raise _ellipticError(Elliptic.fRG, x, y, z, cause=X, txt=t)

927

928 @staticmethod

929 def fRJ(x, y, z, P): # *over

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

931

932 @return: C{RJ(x, y, z, P)}.

933

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

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

936 '''

937 try:

938 return float(_RJ(x, y, z, P))

939 except Exception as X:

940 raise _ellipticError(Elliptic.fRJ, x, y, z, P, cause=X)

941

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

943

944

945class _Ek(Elliptic):

946 '''(INTERNAL) Low-overhead C{Elliptic} for C{Ellipse._Ek}.

947 '''

948 _alpha2 = _0_0

949 _alphap2 = _1_0

950 cE = _0_0 # overide Property_RO

951# cK = _0_0 # ditto

952 _iteration = 0

953

954 def __init__(self, k2):

955 # assert 0 < k2 < 1

956 self._k2 = k2

957 self._kp2 = kp2 = _1_0 - k2

958 self.cE = self._cE(kp2)

959# self.cK = self._cK(kp2)

960

961

962class EllipticError(_ValueError):

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

964 '''

965 pass

966

967

968class Elliptic3Tuple(_NamedTuple):

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

970 '''

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

972 _Units_ = ( Scalar, Scalar, Scalar)

973

974

975class _Hdot(dict):

976 '''(INTERNAL) Caching helper for C{_Horner} and C{_RF3}.

977 '''

978 def __call__(self, F, h, *Xys):

979 k = Xys[1] # unique key

980 ys = self.get(k, None)

981 if ys is None:

982 self[k] = ys = tuple((y / h) for y in Xys[1::2])

983 try:

984 D = Fdot(Xys[0::2], *ys, f2product=True)

985 except (OverflowError, TypeError, ValueError):

986 ts = map(_operator.mul, Xys[0::2], ys)

987 D = Fsum(*ts, nonfinites=True) # _Ksum(0, 0, *ts)

988 return D.fmul(F) # Fsum

989

990_Hdot = _Hdot() # PYCHOK singleton

991

992

993class _List(list):

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

995 '''

996 _a0 = None

997

998 def __init__(self, *xyzp): # x, y, z [, P]

999 list.__init__(self, xyzp)

1000

1001 def a0(self, n):

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

1003 '''

1004 a = list(self)

1005 m = n - len(a)

1006 if m > 0:

1007 a[-1] *= m + 1

1008 self._a0 = a0 = _fsum(a) / n

1009 return a0

1010

1011 def amrs4(self, y, Tol, inst=None):

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

1013 C{lam = fdot(sqrt(x), ... (z))} and C{s = (sqrt(x), ... (P))}.

1014 '''

1015 L = self

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

1017 t = L.threshold(Tol)

1018 m = 1

1019 for i in range(_MAXIT):

1020 d = fabs(a * m)

1021 if d > t: # 3-6 trips

1022 break

1023 s = map2(sqrt, L) # sqrt(x), sqrt(y), sqrt(z) [, sqrt(P)]

1024 Q = _Qdot3(*s) # (s[0] * s[1], s[1] * s[2], s[2] * s[0])

1025 a = Q(a) # An = sum(An, *Q)) / 4

1026 L[:] = map(Q, L) # x = sum(x, *Q) / 4, ...

1027 if y: # yield only if used

1028 r = float(Q) # lam = sum(Q)

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

1030 m *= 4

1031 else: # PYCHOK no cover

1032 raise _ConvergenceError(_MAXIT, d, t, thresh=True)

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

1034 if inst:

1035 inst._iteration += i

1036

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

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

1039 '''

1040 # assert am

1041 a0 = self._a0

1042 _am = _1_over(am)

1043 for x in xs:

1044 yield (a0 - x) * _am

1045

1046 def threshold(self, Tol):

1047 '''Get the convergence C{threshold}.

1048 '''

1049 a0 = self._a0

1050 return max(fabs(x - a0) for x in self) / Tol

1051

1052

1053# class _Qdot3(Fsum):

1054# '''(INTERNAL) "Quarter" 3-dot product.

1055# '''

1056# def __init__(self, x, y, z, *unused): # PYCHOK signature

1057# Fsum.__init__(self, x * y, y * z, z * x)

1058#

1059# def __call__(self, a): # PYCHOK signature

1060# return (self + a).fover(_4_0)

1061

1062

1063class _Qdot3(list):

1064 '''(INTERNAL) "Quarter" 3-dot product.

1065 '''

1066 def __init__(self, x, y, z, *unused): # PYCHOK signature

1067 R = _Rdot(x, y, z, _0_0).partials

1068 list.__init__(self, (0,) + R) # NOT R.fsum2()!

1069

1070 def __call__(self, a):

1071 try:

1072 self[0] = a

1073 q = float(self) * _0_25

1074 finally:

1075 self[0] = 0

1076 return q

1077

1078 def __float__(self):

1079 return _fsum(self) # nonfinites=True

1080

1081

1082def _abm3(x, y, inst=None):

1083 '''(INTERNAL) Yield Carlson 3-tuples C{(xn, yn, m)}.

1084 '''

1085 a, b = sqrt(x), (sqrt(y) if y != _1_0 else y)

1086 if b > a:

1087 b, a = a, b

1088 yield a, -b, _0_5 # (x0 + y0)**2 * _0_5

1089 m = -1

1090 for i in range(_MAXIT):

1091 d = fabs(a - b)

1092 if d <= (a * _TolRG0): # 2..4 trips

1093 break

1094 P = a

1095 a = favg(P, b)

1096 b = sqrt(P * b)

1097 yield a, b, m # (xi - yi)**2 * m

1098 m *= 2

1099 else: # PYCHOK no cover

1100 raise _ConvergenceError(_MAXIT, _over(d, P), _TolRG0)

1101 if inst:

1102 inst._iteration += i

1103 yield a, b, 0 # sentinel: m = 0

1104

1105

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

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

1108 '''

1109 try:

1110 if cn is None or dn is None:

1111 raise ValueError(_invalid_)

1112

1113 if _signBit(cn):

1114 sn, cn = neg_(sn, cn)

1115 r = fX(sn, cn, dn) * _over(PI_2, cX)

1116 return r - atan2(sn, cn)

1117

1118 except Exception as X:

1119 n = NN(_delta_, typename(fX)[1:])

1120 raise _ellipticError(n, sn, cn, dn, cause=X)

1121

1122

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

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

1125 '''

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

1127 return cause, txt, kwds

1128

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

1130

1131 n = typename(where, where)

1132 n = _DOT_(typename(Elliptic), n)

1133 n = _SPACE_(_invokation_, n)

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

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

1136

1137

1138def _Hsum(S, e1, E2, E3, E4, E5, over):

1139 '''(INTERNAL) Horner-like form for C{_RD} and C{_RJ} below.

1140 '''

1141 E22 = E2**2

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

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

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

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

1146 # converted to Horner-like form ...

1147 _1 = _1_0

1148 h = 4084080

1149 S *= e1

1150 S += _Hdot(E5, h, E2, -540540, _1, 471240)

1151 S += _Hdot(E4, h, E2, 612612, E3, -540540, _1, -556920)

1152 S += _Hdot(E3, h, E2, -706860, E22, 675675, E3, 306306, _1, 680680)

1153 S += _Hdot(E2, h, E2, 417690, E22, -255255, _1, -875160)

1154 S += _1

1155 if over:

1156 e1 *= over

1157 return S.fdiv(e1) # Fsum

1158

1159

1160def _3over(a, b):

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

1162 '''

1163 return _over(_3_0, a * b)

1164

1165

1166def _rC(x, y):

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

1168 '''

1169 if x < y: # catch NaN

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

1171 d = y - x

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

1173 elif x == y: # XXX d < EPS0? or EPS02 or _EPSmin

1174 d, r = y, _1_0

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

1176 d = x - y

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

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

1179 d = x - y

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

1181 else: # PYCHOK no cover

1182 d = 0 # y == 0

1183 if d > 0 and x >= 0:

1184 return r / sqrt(d) # float

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

1186

1187

1188def _RD(x, y, z, over=_0_0, inst=None):

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

1190 '''

1191 L = _List(x, y, z)

1192 S = Fsum()

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

1194 if s:

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

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

1197 m *= a

1198 x, y = L.rescale(-m, x, y)

1199 xy = x * y

1200 z = (x + y) / _3_0

1201 z2 = z**2

1202 return _Hsum(S, sqrt(a) * m,

1203 (xy - z2 * _6_0),

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

1205 (xy - z2) * z2 * _3_0,

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

1207

1208

1209def _Rdot(x, y, z, start3):

1210 '''(INTERNAL) "Rotated" C{dot}.

1211 '''

1212 try:

1213 R = Fdot_(x, y, y, z, z, x, start3, _3_0, f2product=True)

1214 except (OverflowError, TypeError, ValueError):

1215 R = Fsum(x * y, y * z, z * x, start3 * _3_0, nonfinites=True)

1216 return R

1217

1218

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

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

1221 '''

1222 for a, b, m in _abm3(x, y, inst): # PYCHOK yield

1223 pass

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

1225

1226

1227def _RF3(x, y, z, inst=None): # 3-arg version

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

1229 '''

1230 L = _List(x, y, z)

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

1232 pass

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

1234 z = neg(x + y)

1235 xy = x * y

1236 e2 = xy - z**2

1237 e3 = xy * z

1238 e4 = e2**2

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

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

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

1242 # converted to Horner-like form ...

1243 _1 = _1_0

1244 h = 240240

1245 S = _Hdot(e3, h, e4, 15015, e3, 6930, e2, -16380, _1, 17160)

1246 S += _Hdot(e2, h, e4, -5775, e2, 10010, _1, -24024)

1247 S += _1

1248 return S.fdiv(sqrt(a)) # Fsum

1249

1250

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

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

1253 '''

1254 rs = [] # len 2..7, incl sentinel

1255 _r = rs.append

1256 for a, b, m in _abm3(x, y, inst): # PYCHOK yield

1257 _r((a - b)**2 * m)

1258 return _over(_fsum(rs) * PI_, a + b) # nonfinites=True, float

1259

1260

1261def _rG3(x, y, z): # 3-arg version # in .triaxials.bases

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

1263 '''

1264 R = _RF3(x, y, z) * z

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

1266 if rd: # Carlson, eq 1.7

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

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

1269 return R.fover(_2_0) # float

1270

1271

1272def _RJ(x, y, z, P, over=_0_0, inst=None):

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

1274 '''

1275 def _xyzp(x, y, z, P):

1276 return (x + P) * (y + P) * (z + P)

1277

1278 L = _List(x, y, z, P)

1279 n = neg(_xyzp(x, y, z, -P))

1280 S = Fsum()

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

1282 if s:

1283 d = _xyzp(*s)

1284 if d:

1285 if n:

1286 r = _rC(_1_0, (n / d**2 + _1_0))

1287 n = n / _64_0 # /= chokes PyChecker

1288 else:

1289 r = _1_0 # == _rC(_1_0, _1_0)

1290 S += r / (d * m)

1291 else: # PYCHOK no cover

1292 return NAN

1293 m *= a

1294 x, y, z = L.rescale(m, x, y, z)

1295 P = neg(Fsum(x, y, z).fover(_2_0))

1296 p2 = P**2

1297 p3 = p2 * P

1298 E2 = _Rdot(x, y, z, -p2)

1299 E2p = E2 * P

1300 xyz = x * y * z

1301 return _Hsum(S.fmul(_6_0), sqrt(a) * m, E2,

1302 Fsum(p3 * _4_0, xyz, E2p, E2p),

1303 Fsum(p3 * _3_0, E2p, xyz, xyz).fmul(P),

1304 p2 * xyz, over) # Fsum

1305

1306

1307class _RJfma(object):

1308 '''(INTERNAL) Carlson, "fma"able.

1309 '''

1310 def __init__(self, *args):

1311 self._Rj = _RJ(*args)

1312

1313 def ma(self, b, c):

1314 r = self._Rj._fma(b, c, nonfinites=True)

1315 # assert r is not self._Rj

1316 return float(r)

1317

1318# **) MIT License

1319#

1321#

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

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

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

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

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

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

1328#

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

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

1331#

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

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

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

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

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

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

1338# OTHER DEALINGS IN THE SOFTWARE.