1

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

3

4 u'''Class L{AuxLat} transcoded to Python from I{Karney}'s C++ class U{AuxLatitude

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AuxLatitude.html>}

6 in I{GeographicLib version 2.2+}.

7

8 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2022-2023) and licensed

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

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

11

12@see: U{Auxiliary latitudes<https:#GeographicLib.SourceForge.io/C++/doc/auxlat.html>}

13 U{On auxiliary latitudes<https:#ArXiv.org/abs/2212.05818>}.

14'''

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

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

17

18 from pygeodesy.auxilats.auxAngle import AuxAngle, AuxBeta, AuxChi, _AuxClass, \

19 AuxMu, AuxPhi, AuxTheta, AuxXi

20 from pygeodesy.auxilats.auxily import Aux, _2cos2x, _sc, _sn, _Ufloats

21 from pygeodesy.basics import isscalar, _reverange, _xinstanceof

22 from pygeodesy.constants import INF, MAX_EXP, MIN_EXP, NAN, NINF, PI_2, PI_4, \

23 _EPSqrt, isfinite, isinf, isnan, _log2, _over, _0_0, \

24 _0_0s, _0_1, _0_25, _0_5, _1_0, _2_0, _3_0, _360_0

25 from pygeodesy.datums import _ellipsoidal_datum, _WGS84, Ellipsoid

26# from pygeodesy.ellipsoids import Ellipsoid # from .datums

27 from pygeodesy.elliptic import Elliptic as _Ef

28 from pygeodesy.errors import AuxError, _xkwds, _xkwds_get, _Xorder

29 from pygeodesy.fmath import cbrt, Fsum

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

31 from pygeodesy.interns import NN, _DOT_, _UNDER_ # _earth_

32 from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS

33 from pygeodesy.props import Property, Property_RO, _update_all

34 from pygeodesy.units import Degrees, Meter

35

36 from math import asinh, atan, atan2, copysign, cosh, fabs, sin, sinh, sqrt

37 try:

38 from math import exp2 as _exp2

39 except ImportError: # Python 3.11-

40

41 def _exp2(x):

42 return pow(_2_0, x)

43

44__all__ = ()

45__version__ = '23.08.05'

46

47_TRIPS = 1024

48

49

50 def _passarg(arg): # PYCHOK to .utily

51 return arg

52

53

54 class AuxLat(AuxAngle):

55 '''Accurate conversion between I{Auxiliary} latitudes on an ellipsoid.

56

57 Latitudes are represented by instances of L{AuxAngle} in order to

58 maintain precision close to the poles, I{Authalic} latitude C{Xi},

59 I{Conformal} C{Chi}, I{Geocentric} C{Theta}, I{Geographic} C{Phi},

60 I{Parametric} C{Beta} and I{Rectifying} C{Mu}.

61

62 @see: I{Karney}'s C++ class U{AuxLatitude

63 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AuxLatitude.html>}.

64 '''

65 _csc = dict() # global coeffs cache: [aL][k], upto max(k) * (4 + 6 + 8) floats

66 _E = _WGS84.ellipsoid

67 _mAL = 6 # 4, 6 or 8 aka Lmax

68

69 def __init__(self, a_earth=_WGS84, f=None, b=None, name=NN, **ALorder):

70 '''New L{AuxLat} instance on an ellipsoid or datum.

71

72 @arg a_earth: Equatorial radius, semi-axis (C{meter}) or an

73 ellipsoid or datum (L{Datum}, L{Ellipsoid},

74 L{Ellipsoid2} or L{a_f2Tuple}).

75 @kwarg f: Flattening: M{(a - b) / a} (C{float}, near zero for

76 spherical), ignored if B{C{a_earth}} is not scalar.

77 @kwarg b: Optional polar radius, semi-axis (C{meter}, same

78 units as B{C{a_earth}}), ignored if B{C{a_earth}}

79 is not scalar.

80 @kwarg ALorder: Optional keyword arguments B{C{ALorder}} to

81 set L{AuxLat}'s C{order}, see property

82 C{ALorder}.

83 @kwarg name: Optional name (C{str}).

84 '''

85 if a_earth is not _WGS84:

86 n = name or AuxLat.__name__

87 try:

88 if b is f is None:

89 E = _ellipsoidal_datum(a_earth, name=n).ellipsoid # XXX raiser=_earth_

90 elif isscalar(a_earth):

91 E = Ellipsoid(a_earth, f=f, b=b, name=_UNDER_(n))

92 else:

93 raise ValueError()

94 except Exception as x:

95 raise AuxError(a_earth=a_earth, f=f, b=b, cause=x)

96 self._E = E

97

98 if name:

99 self.name = name

100 if ALorder:

101 self.ALorder = _xkwds_get(ALorder, ALorder=AuxLat._mAL)

102

103 @Property_RO

104 def a(self):

105 '''Get the equatorial radius (C{Meter}, conventionally).

106 '''

107 return self.ellipsoid.a

108

109 equatoradius = a

110

111 @Property

112 def ALorder(self):

113 '''Get the I{AuxLat} order (C{int}, 4, 6 or 8).

114 '''

115 return self._mAL

116

117 @ALorder.setter # PYCHOK setter!

118 def ALorder(self, order):

119 '''Set the I{AuxLat} order (C{int}, 4, 6 or 8).

120 '''

121 m = _Xorder(_AR2Coeffs, AuxError, ALorder=order)

122 if self._mAL != m:

123 _update_all(self)

124 self._mAL = m

125

126 def _atanhee(self, tphi): # see Ellipsoid._es_atanh, .albers._atanhee

127 # atanh(e * sphi) = asinh(e' * sbeta)

128 f = self.f

129 s = _sn(self._fm1 * tphi) if f > 0 else _sn(tphi)

130 if f: # atanh(e * sphi) = asinh(e' * sbeta)

131 e = self._e

132 s = (atan(e * s) if f < 0 else asinh(self._e1 * s)) / e

133 return s

134

135 def Authalic(self, Phi, **diff_name):

136 '''Convert I{Geographic} to I{Aunthalic} latitude.

137

138 @arg Phi: Geographic latitude (L{AuxAngle}).

139

140 @return: Parametric latitude, C{Beta} (L{AuxAngle}).

141 '''

142 _xinstanceof(AuxAngle, Phi=Phi)

143 # assert Phi._AUX == Aux.PHI

144 tphi = fabs(Phi.tan)

145 if isfinite(tphi) and tphi and self.f:

146 y, x = Phi._yx_normalized

147 q = self._q

148 qv = self._qf(tphi)

149 Dq2 = self._Dq(tphi)

150 Dq2 *= (q + qv) / (fabs(y) + _1_0) # _Dq(-tphi)

151 Xi = AuxXi(copysign(qv, Phi.y), x * sqrt(Dq2),

152 name=_xkwds_get(diff_name, name=Phi.name))

153

154 if _xkwds_get(diff_name, diff=False):

155 if isnan(tphi):

156 d = self._e2m1_sq2

157 else:

158 c = self.Parametric(Phi)._x_normalized

159 d = _over(c, Xi._x_normalized)**3

160 d *= _over(c, x) * _over(_2_0, q)

161 Xi._diff = d

162 else:

163 Xi = AuxXi(*Phi._yx) # diff default

164 # assert Xi._AUX == Aux.XI

165 return Xi

166

167 def AuthalicRadius2(self, exact=False, f_max=_0_1):

168 '''Get the I{Authalic} radius I{squared}.

169

170 @kwarg exact: If C{True}, use the exact expression, otherwise

171 the I{Taylor} series.

172 @kwarg f_max: C{Flattening} not to exceed (C{float}).

173

174 @return: Authalic radius I{squared} (C{meter} I{squared}, same

175 units as the ellipsoid axes).

176

177 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)}

178 exceeds C{f_max}.

179 '''

180 f = self.f

181 if exact or not f:

182 c2 = self.ellipsoid.b2 * self._q # == ellipsoid.c2x * 2

183 elif fabs(f) < f_max:

184 # Using a * (a + b) / 2 as the multiplying factor leads to a rapidly

185 # converging series in n. Of course, using this series isn't really

186 # necessary, since the exact expression is simple to evaluate. However,

187 # we do it for consistency with RectifyingRadius and, presumably, the

188 # roundoff error is smaller compared to that for the exact expression.

189 m = self.ALorder

190 c2 = _MODS.karney._polynomial(self._n, _AR2Coeffs[m], 0, m)

191 c2 *= self.a * (self.a + self.b)

192 else:

193 raise AuxError(exact=exact, f=f, f_max=f_max)

194 return c2 * _0_5

195

196 @Property_RO

197 def b(self):

198 '''Get the polar radius (C{Meter}, conventionally).

199 '''

200 return self.ellipsoid.b # a * (_1_0 - f)

201

202 polaradius = b

203

204 def _coeffs(self, auxout, auxin):

205 # Get the polynomial coefficients as 4-, 6- or 8-tuple

206 aL = self.ALorder # aka Lmax

207 k = Aux._1d(auxout, auxin)

208 try:

209 cs = AuxLat._csc[aL][k]

210 except KeyError:

211 if auxout == auxin:

212 cs = _0_0s(aL) # not cached

213 else:

214 Cx = _CXcoeffs(aL)

215 try:

216 Cx = Cx[auxout][auxin]

217 except KeyError as x:

218 raise AuxError(auxout=auxout, auxin=auxin, cause=x)

219

220 d = x = n = self._n

221 if Aux.use_n2(auxin) and Aux.use_n2(auxout):

222 x = self._n2

223

224 def _m(m):

225 return m // 2

226 else:

227 def _m(m): # PYCHOK expected

228 return m

229

230 i = 0

231 cs = []

232 _p = _MODS.karney._polynomial

233 for r in _reverange(aL):

234 j = i + _m(r) + 1 # order m = j - i - 1

235 cs.append(_p(x, Cx, i, j) * d)

236 d *= n

237 i = j

238 # assert i == len(Cx) and len(cs) == aL

239 AuxLat._csc.setdefault(aL, {})[k] = tuple(cs)

240 return cs

241

242 def Conformal(self, Phi, **diff_name):

243 '''Convert I{Geographic} to I{Conformal} latitude.

244

245 @arg Phi: Geographic latitude (L{AuxAngle}).

246

247 @return: Conformal latitude, C{Chi} (L{AuxAngle}).

248 '''

249 _xinstanceof(AuxAngle, Phi=Phi)

250 # assert Phi._AUX == Aux.PHI

251 tphi = tchi = fabs(Phi.tan)

252 if isfinite(tphi) and tphi and self.f:

253 sig = sinh(self._e2 * self._atanhee(tphi))

254 scsig = _sc(sig)

255 scphi = _sc(tphi)

256 if self.f > 0:

257 # The general expression for tchi is

258 # tphi * scsig - sig * scphi

259 # This involves cancellation if f > 0, so change to

260 # (tphi - sig) * (tphi + sig) / (tphi * scsig + sig * scphi)

261 # To control overflow, write as (sigtphi = sig / tphi)

262 # (tphi - sig) * (1 + sigtphi) / (scsig + sigtphi * scphi)

263 sigtphi = sig / tphi

264 if sig < (tphi * _0_5):

265 tphimsig = tphi - sig

266 else:

267 def _asinh_2(x):

268 return asinh(x) * _0_5

269 # Still possibly dangerous cancellation in tphi - sig.

270 # Write tphi - sig = (1 - e) * Dg(1, e)

271 # Dg(x, y) = (g(x) - g(y)) / (x - y)

272 # g(x) = sinh(x * atanh(sphi * x))

273 # Note sinh(atanh(sphi)) = tphi

274 # Turn the crank on divided differences, substitute

275 # sphi = tphi/_sc(tphi)

276 # atanh(x) = asinh(x/sqrt(1-x^2))

277 e = self._e

278 em1 = self._e2m1 / (_1_0 + e)

279 # assert em1 != 0

280 scb = _sc(self._fm1 * tphi) # sec(beta)

281 scphib = _sc(tphi) / scb # sec(phi) / sec(beta)

282 atphib = _asinh_2(e * tphi / scb) # atanh(e * sphi)

283 atphi = _asinh_2(tphi) # atanh(sphi)

284 t = _asinh_2(em1 * (tphi * scphib)) / em1

285 try:

286 Dg = Fsum(atphi, atphib, t, e * t)

287 except ValueError: # NAN

288 Dg = sum((atphi, atphib, t, e * t))

289 e *= atphib

290 t = atphi - e

291 Dg *= cosh(atphi + e) * _over(sinh(t), t)

292 tphimsig = float(Dg) * em1 # tphi - sig

293 tchi = _over(tphimsig * (_1_0 + sigtphi),

294 scsig + scphi * sigtphi)

295 else:

296 tchi = tphi * scsig - sig * scphi

297

298 n = _xkwds_get(diff_name, name=Phi.name)

299 Chi = AuxChi(tchi, name=n).copyquadrant(Phi)

300

301 if _xkwds_get(diff_name, diff=False):

302 if isinf(tphi): # PYCHOK np cover

303 d = self._conformal_diff

304 else:

305 x = self.Parametric(Phi)._x_normalized

306 d = (self._e2m1 * _over(x, Chi._x_normalized)

307 * _over(x, Phi._x_normalized)) if x else _0_0

308 Chi._diff = d

309 # assrt Chi._AUX == Aux.CHI

310 return Chi

311

312 @Property_RO

313 def _conformal_diff(self):

314 '''(INTERNAL) Constant I{Conformal} diff.

315 '''

316 e = self._e

317 if self.f > 0:

318 ss = sinh(e * asinh(self._e1))

319 d = _over(_1_0, _sc(ss) + ss)

320 else:

321 ss = sinh(-atan(e) * e)

322 d = _sc(ss) - ss

323 return d

324

325 def convert(self, auxout, Zeta_d, exact=False):

326 # Convert I{Auxiliary} or I{scalar} latitude

327 Z = d = Zeta_d

328 if isinstance(Z, AuxAngle):

329 A, auxin = _AuxClass(auxout), Z._AUX

330 if auxin == auxout:

331 pass

332 elif not (isfinite(Z.tan) and Z.tan): # XXX

333 Z = A(*Z._yx, aux=auxout, name=Z.name)

334 elif exact:

335 p = Aux.power(auxout, auxin)

336 if p is None:

337 P = self._fromAux(Z) # Phi

338 Z = self._toAux(auxout, P)

339 Z._iter = P.iteration

340 else:

341 y, x = Z._yx

342 if p:

343 y *= pow(self._fm1, p)

344 Z = A(y, x, aux=auxout, name=Z.name)

345 else:

346 cs = self._coeffs(auxout, auxin)

347 yx = Z._yx_normalized

348 Z = A(*yx, aux=auxout, name=Z.name)

349 # assert Z._yx == yx

350 r = _Clenshaw(True, Z, cs, self.ALorder)

351 Z += AuxAngle.fromRadians(r, aux=auxout)

352 # assert Z._AUX == auxout

353 return Z

354

355 elif isinstance(d, Degrees) or isscalar(d):

356 Z = AuxPhi.fromDegrees(d)

357 d = round((d - Z.toDegrees) / _360_0) * _360_0

358 d += self.convert(auxout, Z, exact).toDegrees

359 return Degrees(d, name=Aux.Greek(auxout))

360

361 raise AuxError(auxout=auxout, Zeta_d=Zeta_d, exact=exact)

362

363 def _Dq(self, tphi):

364 # I{Divided Difference} of (q(1) - q(sphi)) / (1 - sphi).

365 sphi = _sn(tphi)

366 if tphi > 0:

367 scphi = _sc(tphi)

368 d = _1_0 / (scphi**2 * (_1_0 + sphi)) # XXX - sphi

369 if d:

370 # General expression for _Dq(1, sphi) is

371 # atanh(e * d / (1 - e2 * sphi)) / (e * d) +

372 # (1 + e2 * sphi) / ((1 - e2 * sphi * sphi) * e2m1);

373 # atanh( e * d / (1 - e2 * sphi))

374 # = atanh( e * d * scphi/(scphi - e2 * tphi))

375 e2m1, ed = self._e2m1, self._e * d

376 if e2m1 and ed:

377 e2 = self._e2

378 if e2 > 0: # assert self.f > 0

379 scb = _sc(self._fm1 * tphi)

380 q = scphi / scb

381 q *= (scphi + tphi * e2) / (e2m1 * scb)

382 q += asinh(self._e1 * d * scphi / scb) / ed

383 else:

384 s2 = sphi * e2

385 q = (_1_0 + s2) / ((_1_0 - sphi * s2) * e2m1)

386 q += (atan2(ed, _1_0 - s2) / ed) if e2 < 0 else _1_0

387 else:

388 q = INF

389 else:

390 q = self._2_e2m12

391 else: # not reached, open-coded in .Authalic

392 q = _over(self._q - self._qf(tphi), _1_0 - sphi)

393 return q

394

395 @Property_RO

396 def _e(self): # unsigned, (1st) eccentricity

397 return self.ellipsoid.e # sqrt(fabs(self._e2))

398

399 @Property_RO

400 def _e1(self):

401 return sqrt(fabs(self._e12))

402

403 @Property_RO

404 def _e12(self):

405 return _over(self._e2, _1_0 - self._e2)

406

407 @Property_RO

408 def _e12p1(self):

409 return _1_0 / self._e2m1

410

411 @Property_RO

412 def _e2(self): # signed, (1st) eccentricity squared

413 return self.ellipsoid.e2

414

415 @Property_RO

416 def _e2m1(self): # 1 less 1st eccentricity squared

417 return self.ellipsoid.e21 # == ._fm1**2

418

419 @Property_RO

420 def _e2m1_sq2(self):

421 return self._e2m1 * sqrt(self._q * _0_5)

422

423 @Property_RO

424 def _2_e2m12(self):

425 return _2_0 / self._e2m1**2

426

427 @Property_RO

428 def _Ef_fRG_a2b2_PI_4(self):

429 E = self.ellipsoid

430 return _Ef.fRG(E.a2, E.b2) / PI_4

431

432 @Property_RO

433 def ellipsoid(self):

434 '''Get the ellipsoid (L{Ellipsoid}).

435 '''

436 return self._E

437

438 @Property_RO

439 def f(self):

440 '''Get the flattening (C{float}).

441 '''

442 return self.ellipsoid.f

443

444 flattening = f

445

446 @Property_RO

447 def _fm1(self): # 1 - flattening

448 return self.ellipsoid.f1

449

450 def _fromAux(self, Zeta, **name): # no diff

451 '''Convert I{Auxiliary} to I{Geographic} latitude.

452

453 @arg Zeta: Auxiliary latitude (L{AuxAngle}).

454

455 @return: Geographic latitude, I{Phi} (L{AuxAngle}).

456 '''

457 _xinstanceof(AuxAngle, Zeta=Zeta)

458 aux = Zeta._AUX

459 n = _xkwds_get(name, name=Zeta.name)

460 f = self._fromAuxCase.get(aux, None)

461 if f is None:

462 Phi = AuxPhi(NAN, name=n)

463 elif callable(f):

464 t = self._fm1

465 t *= f(t)

466 Phi = _Newton(t, Zeta, self._toZeta(aux), name=n)

467 else: # assert isscalar(f)

468 y, x = Zeta._yx

469 Phi = AuxPhi(y / f, x, name=n)

470 # assert Phi._AUX == Aux.PHI

471 return Phi

472

473 @Property_RO

474 def _fromAuxCase(self):

475 '''(INTERNAL) switch(auxin): ...

476 '''

477 return {Aux.AUTHALIC: cbrt,

478 Aux.CONFORMAL: _passarg,

479 Aux.GEOCENTRIC: self._e2m1,

480 Aux.GEOGRAPHIC: _1_0,

481 Aux.PARAMETRIC: self._fm1,

482 Aux.RECTIFYING: sqrt}

483

484 def Geocentric(self, Phi, **name): # no diff

485 '''Convert I{Geographic} to I{Geocentric} latitude.

486

487 @arg Phi: Geographic latitude (L{AuxAngle}).

488

489 @return: Geocentric latitude, C{Phi} (L{AuxAngle}).

490 '''

491 _xinstanceof(AuxAngle, Phi=Phi)

492 # assert Phi._AUX == Aux.PHI

493 return AuxTheta(Phi.y * self._e2m1, Phi.x,

494 name=_xkwds_get(name, name=Phi.name))

495

496 def Geodetic(self, Phi, **name): # no diff

497 '''Convert I{Geographic} to I{Geodetic} latitude.

498

499 @arg Phi: Geographic latitude (L{AuxAngle}).

500

501 @return: Geodetic latitude, C{Phi} (L{AuxAngle}).

502 '''

503 _xinstanceof(AuxAngle, Phi=Phi)

504 # assert Phi._AUX == Aux.PHI

505 return AuxPhi(Phi, name=_xkwds_get(name, name=Phi.name))

506

507 @Property_RO

508 def _n(self): # 3rd flattening

509 return self.ellipsoid.n

510

511 @Property_RO

512 def _n2(self):

513 return self._n**2

514

515 def Parametric(self, Phi, **name): # no diff

516 '''Convert I{Geographic} to I{Parametric} latitude.

517

518 @arg Phi: Geographic latitude (L{AuxAngle}).

519

520 @return: Parametric latitude, C{Beta} (L{AuxAngle}).

521 '''

522 _xinstanceof(AuxAngle, Phi=Phi)

523 # assert Phi._AUX == Aux.PHI

524 return AuxBeta(Phi.y * self._fm1, Phi.x,

525 name=_xkwds_get(name, name=Phi.name))

526

527 Reduced = Parametric

528

529 @Property_RO

530 def _q(self): # constant _q

531 q, f = self._e12p1, self.f

532 if f:

533 e = self._e

534 q += (asinh(self._e1) if f > 0 else atan(e)) / e

535 else:

536 q += _1_0

537 return q

538

539 def _qf(self, tphi):

540 # function _q: atanh(e * sphi) / e + sphi / (1 - (e * sphi)^2)

541 scb = _sc(self._fm1 * tphi)

542 return self._atanhee(tphi) + (tphi / scb) * (_sc(tphi) / scb)

543

544 def _qIntegrand(self, beta):

545 # pbeta(beta) = integrate(q(beta), beta), with beta in radians

546 # q(beta) = (1-f) * (sin(xi) - sin(chi)) / cos(phi)

547 # = (1-f) * (cos(chi) - cos(xi)) / cos(phi) *

548 # (cos(xi) + cos(chi)) / (sin(xi) + sin(chi))

549 # Fit q(beta)/cos(beta) with Fourier transform

550 # q(beta)/cos(beta) = sum(c[k] * sin((2*k+1)*beta), k, 0, K-1)

551 # then the integral is

552 # pbeta = sum(d[k] * cos((2*k+2)*beta), k, 0, K-1)

553 # where

554 # d[k] = -1/(4*(k+1)) * (c[k] + c[k+1]) for k in 0..K-2

555 # d[K-1] = -1/(4*K) * c[K-1]

556 Beta = AuxBeta.fromRadians(beta)

557 if Beta.x: # and self._fm1:

558 Ax, _cv = Aux, self.convert

559 Phi = _cv(Ax.PHI, Beta, exact=True)

560 schi, cchi = _cv(Ax.CHI, Phi, exact=True)._yx_normalized

561 sxi, cxi = _cv(Ax.XI, Phi, exact=True)._yx_normalized

562 r = (sxi - schi) if fabs(schi) < fabs(cchi) else \

563 _over(_2cos2x(cchi, cxi), (sxi + schi) * _2_0)

564 r *= _over(self._fm1, Phi._x_normalized * Beta._x_normalized)

565 else: # beta == PI_2, PI3_2, ...

566 r = _0_0 # XXX 0 avoids NAN summation exceptions

567 return r

568

569 def Rectifying(self, Phi, **diff_name):

570 '''Convert I{Geographic} to I{Rectifying} latitude.

571

572 @arg Phi: Geographic latitude (L{AuxAngle}).

573

574 @return: Rectifying latitude, C{Mu} (L{AuxAngle}).

575 '''

576 Beta = self.Parametric(Phi)

577 # assert Beta._AUX == Aux.BETA

578 sb, cb = map(fabs, Beta._yx_normalized)

579 a, ka, ka1 = _1_0, self._e2, self._e2m1

580 b, kb, kb1 = self._fm1, -self._e12, self._e12p1

581 if self.f < 0:

582 a, b = b, a

583 ka, kb = kb, ka

584 ka1, kb1 = kb1, ka1

585 sb, cb = cb, sb

586 # now a, b = larger/smaller semiaxis

587 # Beta measured from larger semiaxis

588 # kb, ka = modulus-squared for distance from Beta = 0, pi/2

589 # NB kb <= 0; 0 <= ka <= 1

590 # sa = b*E(Beta, sqrt(kb))

591 # sb = a*E(Beta',sqrt(ka))

592 # 1 - ka * (1 - sb2) = 1 - ka + ka*sb2

593 sb2 = sb**2

594 cb2 = cb**2

595 da2 = ka1 + ka * sb2

596 db2 = _1_0 - kb * sb2

597 # DLMF Eq. 19.25.9

598 my = b * sb * (_Ef.fRF(cb2, db2, _1_0) -

599 kb * sb2 * _Ef.fRD(cb2, db2, _1_0, _3_0))

600 # DLMF Eq. 19.25.10 with complementary angles

601 mx = a * cb * (_Ef.fRF(sb2, da2, _1_0) * ka1 +

602 ka * cb2 * _Ef.fRD(sb2, _1_0, da2, _3_0) * ka1 +

603 ka * sb / sqrt(da2))

604 # my + mx = 2*_Ef.fRG(a*a, b*b) = a*E(e) = b*E(i*e')

605 # mr = a*E(e)*(2/pi) = b*E(i*e')*(2/pi)

606 if self.f < 0:

607 a, b = b, a

608 my, mx = mx, my

609 mr = (my + mx) / PI_2

610 if mr:

611 my = sin(my / mr)

612 mx = sin(mx / mr) # XXX zero?

613 else: # zero Mu

614 my, mx = _0_0, _1_0

615 n = _xkwds_get(diff_name, name=Phi.name)

616 Mu = AuxMu(my, mx, # normalized

617 name=n).copyquadrant(Phi)

618

619 if _xkwds_get(diff_name, diff=False):

620 d, x = _0_0, Beta._x_normalized

621 if x and mr:

622 if Mu.x and Phi.x and not isinf(Phi.tan):

623 d = b / mr * (x / Mu.x)**2 \

624 * (x / Phi._x_normalized)

625 else:

626 d = mr / a

627 Mu._diff = self._fm1 * d

628 return Mu

629

630 def RectifyingRadius(self, exact=False):

631 '''Get the I{Rectifying} radius.

632

633 @arg exact: If C{True}, use the exact expression,

634 otherwise the I{Taylor} series.

635

636 @return: Rectifying radius (L{Meter}, same units

637 as the ellipsoid axes).

638 '''

639 r = self._Ef_fRG_a2b2_PI_4 if exact else self._RectifyingR

640 return Meter(r, name=self.RectifyingRadius.__name__)

641

642 @Property_RO

643 def _RectifyingR(self):

644 m = self.ALorder

645 d = _MODS.karney._polynomial(self._n2, _RRCoeffs[m], 0, m // 2)

646 return d * (self.a + self.b) * _0_5

647

648 def _toAux(self, auxout, Phi, **diff_name):

649 '''Convert I{Geographic} to I{Auxiliary} latitude.

650

651 @arg auxout: I{Auxiliary} kind (C{Aux.KIND}).

652 @arg Phi: Geographic latitude (L{AuxLat}).

653

654 @return: Auxiliary latitude, I{Eta} (L{AuxLat}).

655 '''

656 _xinstanceof(AuxAngle, Phi=Phi)

657 # assert Phi._AUX == Aux.PHI

658 n = _xkwds_get(diff_name, name=Phi.name)

659 m = _toAuxCase.get(auxout, None)

660 if m: # callable

661 A = m(self, Phi, **_xkwds(diff_name, name=n))

662 elif auxout == Aux.GEOGRAPHIC: # == .GEODETIC

663 A = AuxPhi(Phi, name=n)

664 A._diff = _1_0

665 else: # auxout?

666 A = AuxPhi(NAN, name=n)

667 # assert A._AUX == auxout

668 return A

669

670 def _toZeta(self, zetaux):

671 '''Return a (lean) function to create C{AuxPhi(tphi)} and

672 convert that into C{AuxAngle} of (fixed) kind C{zetaux}

673 for use only inside the C{_Newton} loop.

674 '''

675 class _AuxPhy(AuxPhi):

676 # lean C{AuxPhi} instance.

677 # _x = _1_0

678

679 def __init__(self, y): # PYCHOK signature

680 self._y = y

681

682 class _AuxPhy1(_AuxPhy):

683 # lean C{AuxPhi} instance.

684 _diff = _1_0

685

686 m = _toAuxCase.get(zetaux, None)

687 if m: # callable

688

689 def _toZeta(tphi):

690 return m(self, _AuxPhy(tphi), diff=True)

691

692 elif zetaux == Aux.GEOGRAPHIC: # == .GEODETIC

693 _toZeta = _AuxPhy1

694

695 else: # zetaux?

696

697 def _toZeta(unused): # PYCHOK expected

698 return _AuxPhy(NAN)

699

700 return _toZeta

701

702

703# switch(auxout): ...

704_toAuxCase = {Aux.AUTHALIC: AuxLat.Authalic,

705 Aux.CONFORMAL: AuxLat.Conformal,

706 Aux.GEOCENTRIC: AuxLat.Geocentric,

707 Aux.PARAMETRIC: AuxLat.Parametric,

708 Aux.RECTIFYING: AuxLat.Rectifying}

709

710

711 def _Clenshaw(sinp, Zeta, cs, K):

712 sz, cz = Zeta._yx # isnormal

713 # Evaluate sum(c[k] * sin((2*k+2) * Zeta)) if sinp else

714 # sum(c[k] * cos((2*k+2) * Zeta))

715 x = _2cos2x(cz, sz) # 2 * cos(2*Zeta)

716 if isfinite(x):

717 U0, U1 = Fsum(), Fsum()

718 else: # avoid Fsum(NAN) exceptions

719 U0 = U1 = _0_0

720 # assert len(cs) == K

721 for r in _reverange(K):

722 U1 -= U0 * x + cs[r]

723 U1, U0 = U0, -U1

724 # u0*f0(Zeta) - u1*fm1(Zeta)

725 # f0 = sin(2*Zeta) if sinp else cos(2*Zeta)

726 # fm1 = 0 if sinp else 1

727 U0 *= (sz * cz * _2_0) if sinp else (x * _0_5)

728 return float(U0) if sinp else float(U0 - U1)

729

730

731 def _CXcoeffs(aL): # PYCHOK in .auxilats.__main__

732 '''(INTERNAL) Get the C{CX_4}, C{_6} or C{_8} coefficients.

733 '''

734 try: # from pygeodesy.auxilats._CX_x import _coeffs_x as _coeffs

735 _CX_x = _DOT_(_MODS.auxilats.__name__, _UNDER_('_CX', aL))

736 _coeffs = _MODS.getattr(_CX_x, _UNDER_('_coeffs', aL))

737 except (AttributeError, ImportError, KeyError, TypeError) as x:

738 raise AuxError(ALorder=aL, cause=x)

739 # assert _coeffs.ALorder == aL

740 # assert _coeffs.n == Aux.len(aL)

741 return _coeffs

742

743

744 def _Newton(tphi, Zeta, _toZeta, **name):

745 # Newton's method

746 _l2 = _log2

747 tz = fabs(Zeta.tan)

748 ltz = _l2(tz) if tz else NINF

749 if isfinite(ltz):

750 tphi = _over(tz, tphi)

751 ltphi = _l2(tphi)

752 ltmin = min(ltphi, MIN_EXP)

753 ltmax = max(ltphi, MAX_EXP)

754# auxin = Zeta._AUX

755 s, n, __2 = 0, 3, _0_5 # n = i + 2

756 _abs, _ov = fabs, _over

757 _p2, _TOL = _exp2, _EPSqrt

758 for i in range(1, _TRIPS): # 1Ki!

759 # _toAux(auxin, AuxPhi(tphi), diff=True)

760 Z1 = _toZeta(tphi)

761 # assert Z1._AUX == auxin

762 tz1, s_ = Z1.tan, s

763 if tz1 > tz:

764 ltmax, s = ltphi, +1

765 elif tz1 < tz:

766 ltmin, s = ltphi, -1

767 else:

768 break

769 # derivative dtan(Zeta)/dtan(Pphi)

770 # to dlog(tan(Zeta))/dlog(tan(Phi))

771 # Zeta.diff *= tphi/tzeta1

772 d = _ov(tphi, tz1) * Z1.diff

773 d = _ov(ltz - _l2(tz1), d)

774 ltphi += d

775 tphi = _p2(ltphi)

776 if _abs(d) < _TOL:

777 i += 1

778 # _toAux(auxin, AuxPhi(tphi), diff=True)

779 Z = _toZeta(tphi)

780 tphi -= _ov(Z.tan - tz, Z.diff)

781 break

782 if (i > n and (s * s_) < 0) or not ltmin < ltphi < ltmax:

783 s, n = 0, i + 2

784 ltphi = (ltmin + ltmax) * __2

785 tphi = _p2(ltphi)

786 else:

787 i = _TRIPS

788 Phi = AuxPhi(tphi, **name).copyquadrant(Zeta)

789 Phi._iter = i

790 else:

791 Phi = AuxPhi(Zeta, **name) # diff as-as

792 Phi._iter = 0

793 # assert Phi._AUX == Aux.PHI

794 return Phi

795

796

797_f, _u = float, _Ufloats()

798_1__f3 = -1 / _f(3) # XXX +1 / _f(3)

799_AR2Coeffs = {4: _u(4 / _f(315), 4 / _f(105), 4 / _f(15), _1__f3),

800 6: _u(4 / _f(1287), 4 / _f(693), 4 / _f(315), 4 / _f(105),

801 4 / _f(15), _1__f3),

802 8: _u(4 / _f(3315), 4 / _f(2145), 4 / _f(1287), 4 / _f(693),

803 4 / _f(315), 4 / _f(105), 4 / _f(15), _1__f3)}

804_RRCoeffs = {4: _u(1 / _f(64), _0_25),

805 6: _u(1 / _f(256), 1 / _f(64), _0_25),

806 8: _u(25 / _f(16384), 1 / _f(256), 1 / _f(64), _0_25)} # PYCHOK used!

807 del _f, _u, _Ufloats, _1__f3

808# assert set(_AR2Coeffs.keys()) == set(_RRCoeffs.keys())

809

810__all__ += _ALL_DOCS(AuxLat)

811

812# **) MIT License

813#

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

815#

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

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

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

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

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

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

822#

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

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

825#

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

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

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

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

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

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

832# OTHER DEALINGS IN THE SOFTWARE.