Coverage for pygeodesy/elliptic.py: 96%

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 the main 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 by 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, _EPSjam, _EPStol, \

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

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

73from pygeodesy.errors import _ConvergenceError, _ValueError

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

75from pygeodesy.fsums import Fsum, _fsum

76from pygeodesy.internals import _Enum, typename

77from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _invalid_, _negative_

78from pygeodesy.karney import _K_2_0, _norm180, _remainder2, _signBit, _sincos2

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

80from pygeodesy.named import _Named, _NamedTuple, unstr

81from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all

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

83from pygeodesy.units import Scalar, Scalar_

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

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

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

88# import operator as _operator # from .fmath

90__all__ = _ALL_LAZY.elliptic

91__version__ = '26.11.09'

93_TolRD = zqrt(EPS * 0.002)

94_TolRF = zqrt(EPS * 0.030)

95_TolRG0 = _EPStol * 2.7

96_MAXIT = 32 # 6..18 sufficient

99class _Cs(_Enum):

100 '''(INTERAL) Complete Integrals cache.

101 '''

102 pass # PYLINT W0107

103

104

105class Elliptic(_Named):

106 '''Elliptic integrals and functions.

107

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

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

110 '''

111# _alpha2 = 0

112# _alphap2 = 0

113# _eps = EPS

114# _k2 = 0

115# _kp2 = 0

116

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

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

119

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

121

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

123

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

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

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

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

128 '''

129 if name:

130 self.name = name

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

132

133 @Property_RO

134 def alpha2(self):

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

136 '''

137 return self._alpha2

138

139 @Property_RO

140 def alphap2(self):

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

142 '''

143 return self._alphap2

144

145 @Property_RO

146 def _Bcd2mn4(self):

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

148 '''

149 a = _1_0

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

151 if _signBit(b): # PYCHOK no cover

152 d = a - b

153 b = neg(b / d)

154 d = sqrt(d)

155 c, mn = _sentuple2(_Bmns(a, b))

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

157 self._iteration += len(mn)

158 return c, d, cd, mn

159

160 def _Bsncndn3(self, x):

161 '''(INTERNAL) Bulirsch' C{sncndn} routine.

162 '''

163 if self.kp2:

164 # Numerische Mathematik 7, 1965, P 89

165 # implements DLMF Eqs 22.17.2 - 22.17.4

166 c, d, cd, mn = self._Bcd2mn4

167 sn, cn = _sincos2(cd * x)

168 dn = _1_0

169 if sn:

170 a = cn / sn

171 c *= a

172 for m, n in mn:

173 a *= c

174 c *= dn

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

176 a = c / m

177 a = _1_0 / hypot1(c)

178 sn = _flipsign(a, sn)

179 cn = sn * c

180 if d: # and _signBit(self.kp2): # PYCHOK no cover

181 cn, dn = dn, cn

182 sn = sn / d # /= chokes PyChecker

183 else: # .kp2 == 0

184 sn, cn, dn = _sncndn3(x)

185 return sn, cn, dn

186

187 @Property_RO

188 def cD(self):

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

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

191 '''

192 return self._cDEKEeps.cD # PYLINT E1101

193

194 @Property_RO

195 def _cDEKEeps(self):

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

197 '''

198 k2, kp2 = self.k2, self.kp2

199 if k2:

200 if kp2:

201 try: # self._iteration = 0

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

203 cE = self._cE(kp2)

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

205 cK = self._cK(kp2)

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

207 # above at <https://DLMF.NIST.gov/19.25.E1>

208 cD = _rD_3(kp2, _1_0)

209 cKE = k2 * cD

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

211

212 except Exception as X:

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

214 else:

215 cD = cK = cKE = INF

216 cE = _1_0

217 eps = k2

218 else:

219 cD = PI_4

220 cE = cK = PI_2

221 cKE = _0_0 # k2 * cD

222 eps = EPS

223

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

225

226 @Property_RO

227 def cE(self):

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

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

230 '''

231 return self._cDEKEeps.cE # PYLINT E1101

232

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

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

235

236 @Property_RO

237 def cG(self):

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

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

240 '''

241 return self._cGHPi.cG # PYLINT E1101

242

243 @Property_RO

244 def _cGHPi(self):

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

246 '''

247 alpha2, alphap2, kp2 = self.alpha2, self.alphap2, self.kp2

248 try:

249 # self._iteration = 0

250 if alpha2:

251 if alphap2:

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

253 cK = self.cK

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

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

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

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

258 else:

259 cG = cH = _rC(_1_0, alphap2)

260 cPi = INF # XXX or NAN?

261 else:

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

263 else:

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

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

266 # So write (for alpha2 = 0)

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

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

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

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

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

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

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

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

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

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

277 cH = PI_4 if kp2 == 1 else (

278 _1_0 if kp2 == 0 else (_rD_3(_1_0, kp2) * kp2))

279

280 except Exception as X:

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

282 alphap2=alphap2, cause=X)

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

284

285 @Property_RO

286 def cH(self):

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

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

289 '''

290 return self._cGHPi.cH # PYLINT E1101

291

292 @Property_RO

293 def cK(self):

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

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

296 '''

297 return self._cDEKEeps.cK # PYLINT E1101

298

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

300 return _rF2(kp2, _1_0, self)

301

302 @Property_RO

303 def cKE(self):

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

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

306 '''

307 return self._cDEKEeps.cKE # PYLINT E1101

308

309 @Property_RO

310 def cPi(self):

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

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

313 '''

314 return self._cGHPi.cPi # PYLINT E1101

315

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

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

318

319 @arg sn: sin(φ).

320 @arg cn: cos(φ).

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

322

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

324

325 @raise EllipticError: Invalid invokation or no convergence.

326 '''

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

328

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

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

331

332 @arg sn: sin(φ).

333 @arg cn: cos(φ).

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

335

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

337

338 @raise EllipticError: Invalid invokation or no convergence.

339 '''

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

341

342 def deltaEinv(self, stau, ctau):

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

344

345 @arg stau: sin(τ)

346 @arg ctau: cos(τ)

347

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

349

350 @raise EllipticError: No convergence.

351 '''

352 try:

353 if _signBit(ctau): # pi periodic

354 stau, ctau = neg_(stau, ctau)

355 t = atan2(stau, ctau)

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

357

358 except Exception as X:

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

360

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

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

363

364 @arg sn: sin(φ).

365 @arg cn: cos(φ).

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

367

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

369

370 @raise EllipticError: Invalid invokation or no convergence.

371 '''

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

373

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

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

376

377 @arg sn: sin(φ).

378 @arg cn: cos(φ).

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

380

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

382

383 @raise EllipticError: Invalid invokation or no convergence.

384 '''

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

386

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

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

389

390 @arg sn: sin(φ).

391 @arg cn: cos(φ).

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

393

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

395

396 @raise EllipticError: Invalid invokation or no convergence.

397 '''

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

399

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

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

402

403 @arg sn: sin(φ).

404 @arg cn: cos(φ).

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

406

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

408 (C{float}).

409

410 @raise EllipticError: Invalid invokation or no convergence.

411 '''

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

413

414 def _Einv(self, x):

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

416 '''

417 E2 = self.cE * _2_0

418 n = floor(x / E2 + _0_5)

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

420 # linear approximation

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

422 Phi = Fsum(phi)

423 # first order correction

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

425 # self._iteration = 0

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

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

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

429 _Phi2 = Phi.fsum2f_ # aggregate

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

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

432 if dn:

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

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

435 else: # PYCHOK no cover

436 d = _0_0 # XXX continue?

437 if fabs(d) < _EPStol: # 3-4 trips

438 self._iteration += i

439 break

440 else: # PYCHOK no cover

441 raise _ConvergenceError(_MAXIT, d, _EPStol)

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

443

444 @Property_RO

445 def eps(self):

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

447 '''

448 return self._cDEKEeps.eps # PYLINT E1101

449

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

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

452 Jacobi elliptic functions.

453

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

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

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

457

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

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

460

461 @raise EllipticError: Invalid invokation or no convergence.

462 '''

463 def _fD(sn, cn, dn):

464 r = fabs(sn)**3

465 if r:

466 R = _RD(cn**2, dn**2, _1_0, _3_0 / r, self)

467 r = float(R)

468 return r

469

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

471 self.deltaD, _fD)

472

473 def fDelta(self, sn, cn):

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

475

476 @arg sn: sin(φ).

477 @arg cn: cos(φ).

478

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

480 '''

481 try:

482 k2 = self.k2

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

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

485 return sqrt(s) if s else _0_0

486

487 except Exception as X:

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

489

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

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

492 Jacobi elliptic functions.

493

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

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

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

497

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

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

500

501 @raise EllipticError: Invalid invokation or no convergence.

502 '''

503 def _fE(sn, cn, dn):

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

505 '''

506 if sn:

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

508 kp2, k2 = self.kp2, self.k2

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

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

511 if k2:

512 Ei -= _RD(cn2, dn2, _1_0, _3_over(k2, sn**2), self)

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

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

515 if kp2:

516 Ei += (_RD( cn2, _1_0, dn2, _3_over(k2, sn**2), self) +

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

518 else: # PYCHOK no cover

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

520 Ei -= _RD(dn2, _1_0, cn2, _3_over(kp2, sn**2), self)

521 Ei *= fabs(sn)

522 else: # PYCHOK no cover

523 Ei = _0_0

524 return float(Ei)

525

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

527 self.deltaE, _fE,

528 k2_0=self.k2==0)

529

530 def fEd(self, deg):

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

532 the argument given in C{degrees}.

533

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

535

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

537

538 @raise EllipticError: No convergence.

539 '''

540 if _K_2_0:

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

542 elif fabs(deg) < _180_0:

543 e = _0_0

544 else:

545 e = ceil(deg / _360_0 - _0_5)

546 deg -= e * _360_0

547 e *= self.cE * _4_0

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

549

550 def fEinv(self, x):

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

552

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

554

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

556 (C{float}).

557

558 @raise EllipticError: No convergence.

559 '''

560 try:

561 return self._Einv(x)

562 except Exception as X:

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

564

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

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

567 Jacobi elliptic functions.

568

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

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

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

572

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

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

575

576 @raise EllipticError: Invalid invokation or no convergence.

577 '''

578 def _fF(sn, cn, dn):

579 r = fabs(sn)

580 if r:

581 R = _RF3(cn**2, dn**2, _1_0, self)

582 r = float(R.fmul(r))

583 return r

584

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

586 self.deltaF, _fF,

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

588

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

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

591 Jacobi elliptic functions.

592

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

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

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

596

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

598

599 @raise EllipticError: Invalid invokation or no convergence.

600

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

602 geodesic in terms of this combination of elliptic

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

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

605 riIOAAAAQAAJ&pg=PA181>}.

606

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

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

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

610 '''

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

612 self.cG, self.deltaG)

613

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

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

616 Jacobi elliptic functions.

617

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

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

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

621

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

623

624 @raise EllipticError: Invalid invokation or no convergence.

625

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

627 on the geodesic in terms of this combination of

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

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

630

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

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

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

634 '''

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

636 self.cH, self.deltaH)

637

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

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

640 Jacobi elliptic functions.

641

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

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

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

645

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

647

648 @raise EllipticError: Invalid invokation or no convergence.

649 '''

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

651 dn = self.fDelta(phi_or_sn, cn) # in .triaxials.triaxial5

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

653 self.cPi, self.deltaPi)

654

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

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

657 '''

658 def _fX(sn, cn, dn):

659 if sn:

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

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

662 if aX:

663 sn2 = sn**2

664 P = sn2 * self.alphap2 + cn2

665 R += _RJ(cn2, dn2, _1_0, P, _3_over(aX, sn2), self)

666 R *= fabs(sn)

667 else: # PYCHOK no cover

668 R = _0_0

669 return float(R)

670

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

672

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

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

675 '''

676 # self._iteration = 0 # aggregate

677 phi = sn = phi_or_sn

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

679 if k2_0: # C++ version 2.4

680 return phi

681 elif kp2_0:

682 return asinh(tan(phi))

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

684 if fabs(phi) >= PI:

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

686 # fall through

687 elif cn is None or dn is None:

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

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

690

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

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

693 else:

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

695 return copysign0(xi, sn)

696

697# def _Jam(self, x):

698# '''Jacobi amplitude function.

699#

700# @return: C{phi} (C{float}).

701# '''

702# _, _, _, phi = self._Jsncndn4(x)

703# return phi

704

705 @Property_RO

706 def _Jpac2(self):

707 '''(INTERNAL) Get Jacobi's 2-tuple C{(p, ac)}.

708 '''

709 p = a = _1_0

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

711 # assert b and c

712 if c < 0: # Sala Eq. 5.8

713 p = sqrt(b)

714 b = _1_over(b)

715# c *= b # unused

716 a, ac = _sentuple2(_Jacs(a, b, c))

717 i = len(ac)

718 p *= 2**i * a

719 self._iteration += i

720 return p, ac

721

722 def _Jsncndn4(self, x):

723 '''(INTERNAL) Jacobi's C{(sn, cn, dn, phi)}, C++ v 2.8.

724 '''

725 if self.k2:

726 if self.kp2: # reduce C{x} as C{y}

727 y, n = _remainder2(x, self.cK)

728 # now x = n*2 + cK + y where y in [-cK, cK], thus

729 # am(x) = am(y) + n*PI with am(y) in [-PI_2, PI_2]

730 if not y:

731 sn = phi = y

732 cn = dn = _1_0

733 elif fabs(y) == self.cK:

734 cn, dn = _0_0, sqrt(self.kp2)

735 sn, phi = _1_0, PI_2

736 if y < 0:

737 sn, phi = (-sn), -phi

738 else: # see DLMF Sec 22.20(ii) and U{Sala 1989

739 # <https://doi.org/10.1137/0520100>}, Sec 5.

740 p, ac = self._Jpac2

741 phi = _Jphi(p * y, ac, self.k2)

742 # use cn = 0 if fabs(phi) > PI_2

743 sn, cn, dn = self._sncndn3(phi, True)

744 if (n % 2): # isodd(n)

745 cn, sn = (-cn), -sn

746 phi += PI * n

747 else: # .kp2 == 0

748 sn, cn, dn = _sncndn3(x)

749 phi = atan(sinh(x)) # gd(x)

750 else: # .k2 == 0

751 sn, cn = _sincos2(x)

752 dn, phi = _1_0, x

753 return sn, cn, dn, phi # XXX phi unused

754

755 @Property_RO

756 def k2(self):

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

758 '''

759 return self._k2

760

761 @Property_RO

762 def kp2(self):

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

764 '''

765 return self._kp2

766

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

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

769

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

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

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

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

774

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

776 or B{C{alphap2}}.

777

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

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

780 that these conditions are met to enable accuracy to be

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

782 '''

783 _update_all(self, Base=Property_RO, needed=5)

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

785

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

787 '''(INITERNAL) Reset this elliptic.

788 '''

789 def _1p2(kp2, k2):

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

791

792 def _S(**kwds):

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

794

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

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

797

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

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

800

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

802 # K E D

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

804 # k = 1: inf 1 inf

805 # Pi G H

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

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

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

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

810 #

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

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

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

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

815 # Pi(alpha2, 1) = inf

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

817

818 self._iteration = 0

819

820 def sncndn(self, x, jam=False): # XXX jam=True tests

821 '''The Jacobi amplitude or Bulirsch' elliptic function.

822

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

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

825 Bulirsch' elliptic function (C{bool}).

826

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

828

829 @raise EllipticError: No convergence.

830 '''

831 i = self._iteration

832 try: # Jacobi amplitude, C++ am v 2.8

833 sn_cn_dn = self._Jsncndn4(x)[:3] if jam else \

834 self._Bsncndn3(x) # handles .kp2 == 0

835 except Exception as X:

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

837

838 i = self._iteration - i

839 return Elliptic3Tuple(*sn_cn_dn, iteration=i or None)

840

841 def _sncndn3(self, phi, cn0=False):

842 '''(INTERNAL) Helper for C{.fEinv}, C{._fXf} and C{._Jsncndn4}.

843 '''

844 sn, cn = _sincos2(phi)

845 if cn0 and _signBit(cn):

846 cn = _0_0

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

848

849 @staticmethod

850 def fRC(x, y):

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

852

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

854

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

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

857 '''

858 return _rC(x, y)

859

860 @staticmethod

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

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

863

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

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

866

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

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

869 '''

870 try:

871 return _rD( y, z, *over) if x == 0 else float(

872 _RD(x, y, z, *over))

873 except Exception as X:

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

875

876 @staticmethod

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

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

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

880

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

882

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

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

885 '''

886 try:

887 return _rF2(y, z) if x == 0 else (

888 _rF2(z, x) if y == 0 else (

889 _rF2(x, y) if z == 0 else (

890 _rC( x, y) if y == z else (

891 _rC( y, z) if z == x else (

892 _rC( z, x) if x == y else float(_RF3(x, y, z)))))))

893 except Exception as X:

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

895

896 @staticmethod

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

898 y = _1_0 - m

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

900 R = _RF3(x, y, z)

901 if m:

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

903 except Exception as X:

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

905 return float(R)

906

907 @staticmethod

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

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

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

911

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

913

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

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

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

917 EllipticFunction_8cpp_source.html#l00096>} version 2.3.

918 '''

919 try:

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

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

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

923 except Exception as X:

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

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

926

927 @staticmethod

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

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

930

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

932

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

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

935 '''

936 try:

937 return float(_RD(y, z, P) if x == P else (

938 _RD(z, x, P) if y == P else (

939 _RD(x, y, P) if z == P else _RJ(x, y, z, P))))

940 except Exception as X:

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

942

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

944

945

946class _Ek(Elliptic):

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

948 '''

949 _alpha2 = _0_0

950 _alphap2 = _1_0

951 cE = _0_0 # overide Property_RO

952# cK = _0_0 # ditto

953 _iteration = 0

954

955 def __init__(self, k2): # PYLINT W0231

956 # assert 0 < k2 < 1

957 self._k2 = k2

958 self._kp2 = kp2 = _1_0 - k2

959 self.cE = self._cE(kp2)

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

961

962

963class EllipticError(_ValueError):

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

965 '''

966 pass # PYLINT W0107

967

968

969class Elliptic3Tuple(_NamedTuple):

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

971 '''

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

973 _Units_ = ( Scalar, Scalar, Scalar)

974

975

976class _Hdot(dict):

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

978 '''

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

980 k = Xys[1] # unique key

981 ys = self.get(k, None)

982 if ys is None:

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

984 try:

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

986 except (OverflowError, TypeError, ValueError):

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

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

989 return D.fmul(F) # Fsum

990

991_Hdot = _Hdot() # PYCHOK singleton

992

993

994class _List(list):

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

996 '''

997 _a0 = None

998

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

1000 list.__init__(self, xyzp)

1001

1002 def a0(self, n):

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

1004 '''

1005 a = list(self)

1006 m = n - len(a)

1007 if m > 0:

1008 a[-1] *= m + 1

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

1010 return a0

1011

1012 def amrs4(self, yd, Tol, inst=None):

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

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

1015 '''

1016 L = self

1017 a = L.a0(5 if yd else 3) # A0

1018 t = L.threshold(Tol)

1019 m = 1

1020 for i in range(_MAXIT):

1021 d = fabs(a * m)

1022 if d > t: # 3-6 trips

1023 break

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

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

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

1027 r = L[2] # "previous" z, s[2] == sqrt(z)

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

1029 if yd: # yield if used

1030 r += float(Q) # lam = sum(*Q)

1031 yield a, m, r, s

1032 m *= 4

1033 else: # PYCHOK no cover

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

1035 if inst:

1036 inst._iteration += i

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

1038

1039 def rescale(self, am, *xyz):

1040 '''Rescale C{x}, C{y}, ... C{z}.

1041 '''

1042 # assert am

1043 a0 = self._a0

1044 return (_over(a0 - x, am) for x in xyz)

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 t = _Rdot(x, y, z, _0_0).partials

1068 list.__init__(self, (0,) + t) # 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

1082class _RJfma(object):

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

1084 '''

1085 def __init__(self, *args):

1086 self._Rj = _RJ(*args)

1087

1088 def ma(self, b, c):

1089 '''Return C{float(._Rj * b + c)}.

1090 '''

1091 R = self._Rj._fma(b, c, nonfinites=True)

1092 # assert R is not self._Rj

1093 return float(R)

1094

1095

1096def _abm3(x, y, inst=None, Tol=_TolRG0):

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

1098 '''

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

1100 if b > a:

1101 b, a = a, b

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

1103 m = -1

1104 for i in range(_MAXIT):

1105 d = fabs(a - b)

1106 if d <= (a * Tol): # 2..4 trips

1107 break

1108 a_ = a

1109 a = favg(a_, b)

1110 b = sqrt(a_ * b)

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

1112 m *= 2

1113 else: # PYCHOK no cover

1114 raise _ConvergenceError(_MAXIT, _over(d, a_), Tol)

1115 if inst:

1116 inst._iteration += i

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

1118

1119

1120def _Bmns(a, b, Tol=_EPStol):

1121 '''(INTERNAL) Yield Bulirsch' C{mn} pairs plus sentinel C{c}.

1122 '''

1123 a_ = a

1124 for _ in range(_MAXIT): # GEOGRAPHICLIB_PANIC

1125 b = sqrt(a_ * b)

1126 yield (a, b)

1127 c = favg(a, b)

1128 d = fabs(a - b)

1129 if d <= (a * Tol): # 3..6 trips, quadratic

1130 break

1131 a_, a = a, c

1132 else: # PYCHOK no cover

1133 raise _ConvergenceError(_MAXIT, _over(d, a_), Tol)

1134 yield c # sentinel

1135

1136

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

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

1139 '''

1140 try:

1141 if cn is None or dn is None:

1142 raise ValueError(_invalid_)

1143

1144 if _signBit(cn):

1145 sn, cn = neg_(sn, cn)

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

1147 return r - atan2(sn, cn)

1148

1149 except Exception as X:

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

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

1152

1153

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

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

1156 '''

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

1158 return cause, txt, kwds

1159

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

1161

1162 n = typename(where, where)

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

1164# n = _SPACE_('invokation', n)

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

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

1167

1168

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

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

1171 '''

1172 E22 = E2**2

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

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

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

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

1177 # converted to Horner-like form ...

1178 _1 = _1_0

1179 h = 4084080

1180 S *= e1

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

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

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

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

1185 S += _1

1186 if over:

1187 e1 *= over

1188 return S.fdiv(e1) # Fsum

1189

1190

1191def _Jacs(a, b, c, Tol=_EPSjam):

1192 '''(INTERNAL) Yield Jacobi's C{ac} pairs plus sentinel C{a}.

1193 '''

1194 # yield (a, sqrt(c)) # unused

1195 for _ in range(_MAXIT): # GEOGRAPHICLIB_PANIC

1196 b = sqrt(a * b)

1197 c = favg(a, -b)

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

1199 # assert a and c

1200 yield (a, c)

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

1202 break

1203 else: # PYCHOK no cover

1204 raise _ConvergenceError(_MAXIT, _over(c, a), Tol)

1205 yield a # sentinel

1206

1207

1208def _Jphi(p, ac, k2):

1209 '''(INTERNAL) Return Jacobi's C{phi}.

1210 '''

1211 for a, c in ac:

1212 p_ = p

1213 a = sin(p) * c / a

1214 p = favg(p, asin(a))

1215 if k2 < 0: # Sala Eq. 5.8

1216 p = p_ - p

1217 return p

1218

1219

1220def _3_over(a, b):

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

1222 '''

1223 return _over(_3_0, a * b)

1224

1225

1226def _rC(x, y):

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

1228 '''

1229 if x < y: # catch NaN

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

1231 d = y - x

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

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

1234 d, r = y, _1_0

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

1236 d = x - y

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

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

1239 d = x - y

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

1241 else: # PYCHOK no cover

1242 d = 0 # y == 0

1243 if d > 0 and x >= 0:

1244 return r / sqrt(d) # PYLINT E0606 float

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

1246

1247

1248def _rD(y, z, *over):

1249 '''(INTERNAL) _RD(x, y, z) for C{x} == 0.

1250 '''

1251 r = _rD_3(y, z) * _3_0

1252 return _over(r, *over) if over else r

1253

1254

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

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

1257 '''

1258 L = _List(x, y, z)

1259 S = Fsum()

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

1261 if s: # except sentinel

1262 S += _3_over(r, s[2] * m) # s[2] = sqrt(z)

1263 m *= a # PYLINT W0631

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

1265 xy = x * y

1266 z = (x + y) / _3_0

1267 z2 = z**2

1268 return _Hsum(S, sqrt(a) * m, # PYLINT W0631

1269 (xy - z2 * _6_0),

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

1271 (xy - z2) * z2 * _3_0,

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

1273

1274

1275def _rD_3(y, z):

1276 '''(INTERNAL) _RD(x, y, z, over=3) for C{x} == 0.

1277 '''

1278 r = y - z

1279 return Fdot_(_rG2(y, z), _2_0,

1280 _rF2(y, z), -z).fover(r * z) if r else \

1281 _over(PI_4, sqrt3(y))

1282

1283

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

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

1286 '''

1287 try:

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

1289 except (OverflowError, TypeError, ValueError):

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

1291 return R # Fsum

1292

1293

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

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

1296 '''

1297 for a, b, _ in _abm3(x, y, inst): # PYCHOK yield

1298 pass

1299 return _over(PI, a + b) # PYLINT W0631 float

1300

1301

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

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

1304 '''

1305 L = _List(x, y, z)

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

1307 pass # only sentinel

1308 x, y = L.rescale(a * m, x, y) # PYLINT W0631

1309 z = neg(x + y)

1310 xy = x * y

1311 e2 = xy - z**2

1312 e3 = xy * z

1313 e4 = e2**2

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

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

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

1317 # converted to Horner-like form ...

1318 _1 = _1_0

1319 h = 240240

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

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

1322 S += _1

1323 return S.fdiv(sqrt(a)) # PYLINT W0631 Fsum

1324

1325

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

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

1328 '''

1329 a = b = None # silence PYLINT

1330 rs = []

1331 _r = rs.append

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

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

1334 rs = _fsum(rs)

1335 return _over(rs * PI_, a + b)

1336

1337

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

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

1340 '''

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

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

1343 if rd: # Carlson, eq 1.7

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

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

1346 return R.fover(_2_0) # float

1347

1348

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

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

1351 '''

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

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

1354

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

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

1357 S = Fsum()

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

1359 if s: # except sentinel

1360 d = _xyzp(*s)

1361 if d:

1362 if n:

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

1364 n = n / _64_0 # /= chokes PyChecker

1365 else:

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

1367 S += r / (d * m)

1368 else: # PYCHOK no cover

1369 return NAN

1370 m *= a # PYLINT W0631

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

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

1373 p2 = P**2

1374 p3 = p2 * P

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

1376 E2p = E2 * P

1377 xyz = x * y * z

1378 return _Hsum(S.fmul(_6_0), sqrt(a) * m, E2, # PYLINT W0631

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

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

1381 p2 * xyz, over) # Fsum

1382

1383

1384def _sentuple2(pairs_sentinel):

1385 '''(INTERNAL) Split sentinel and tuple, I{reversed}.

1386 '''

1387 t = list(pairs_sentinel)

1388 s = t.pop()

1389 # assert isinstance(s, float)

1390 return s, tuple(reversed(t))

1391

1392

1393def _sncndn3(x):

1394 '''(INTERNAL) Return 3-tuple C{(sn, cn, dn)}.

1395 '''

1396 s = tanh(x) # accurate for large abs(x)

1397 c = _1_over(cosh(x))

1398 return s, c, c

1399

1400# **) MIT License

1401#

1403#

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

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

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

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

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

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

1410#

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

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

1413#

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

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

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

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

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

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

1420# OTHER DEALINGS IN THE SOFTWARE.