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:Karney@Alum.MIT.edu>} (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, _sc, _sn, _Ufloats

21 from pygeodesy.basics import isscalar, _reverange, _xinstanceof

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

23 _0_0, _0_0s, _0_1, _0_25, _0_5, _1_0, _2_0, _3_0, \

24 _360_0, isfinite, isinf, isnan, _log2, _over

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 # from .karney

30 from pygeodesy.fsums import Fsum, _sum

31 from pygeodesy.karney import _2cos2x, _polynomial, _ALL_DOCS, cbrt, _MODS

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

33# from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS # from .karney

34 from pygeodesy.props import Property, Property_RO, _update_all

35 from pygeodesy.units import Degrees, Meter

36# from pygeodesy.utily import _passarg # _MODS

37

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

39 try:

40 from math import exp2 as _exp2

41 except ImportError: # Python 3.11-

42

43 def _exp2(x):

44 return pow(_2_0, x)

45

46__all__ = ()

47__version__ = '23.09.18'

48

49_TRIPS = 1024 # XXX 2 or 3?

50

51

52 class AuxLat(AuxAngle):

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

54

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

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

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

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

59

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

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

62 '''

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

64 _E = _WGS84.ellipsoid

65# _Lmax = 0 # overwritten below

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

67

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

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

70

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

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

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

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

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

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

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

78 is not scalar.

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

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

81 C{ALorder}.

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

83 '''

84 if a_earth is not _WGS84:

85 n = name or AuxLat.__name__

86 try:

87 if b is f is None:

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

89 elif isscalar(a_earth):

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

91 else:

92 raise ValueError()

93 except Exception as x:

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

95 self._E = E

96

97 if name:

98 self.name = name

99 if ALorder:

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

101

102 @Property_RO

103 def a(self):

104 '''Get the C{ellipsoid}'s equatorial radius (C{meter}, conventionally).

105 '''

106 return self.ellipsoid.a

107

108 equatoradius = a

109

110 @Property

111 def ALorder(self):

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

113 '''

114 return self._mAL

115

116 @ALorder.setter # PYCHOK setter!

117 def ALorder(self, order):

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

119 '''

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

121 if self._mAL != m:

122 _update_all(self)

123 self._mAL = m

124

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

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

127 f = self.f

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

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

130 e = self._e

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

132 return s

133

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

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

136

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

138

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

140 '''

141 _xinstanceof(AuxAngle, Phi=Phi)

142 # assert Phi._AUX == Aux.PHI

143 tphi = fabs(Phi.tan)

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

145 y, x = Phi._yx_normalized

146 q = self._q

147 qv = self._qf(tphi)

148 Dq2 = self._Dq(tphi)

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

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

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

152

153 if _xkwds_get(diff_name, diff=False):

154 if isnan(tphi):

155 d = self._e2m1_sq2

156 else:

157 c = self.Parametric(Phi)._x_normalized

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

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

160 Xi._diff = d

161 else:

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

163 # assert Xi._AUX == Aux.XI

164 return Xi

165

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

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

168

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

170 the I{Taylor} series.

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

172

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

174 units as the ellipsoid axes).

175

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

177 exceeds C{f_max}.

178 '''

179 f = self.f

180 if exact or not f:

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

182 elif fabs(f) < f_max:

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

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

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

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

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

188 m = self.ALorder

189 c2 = _polynomial(self._n, _AR2Coeffs[m], 0, m)

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

191 else:

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

193 return c2 * _0_5

194

195 @Property_RO

196 def b(self):

197 '''Get the C{ellipsoid}'s polar radius (C{meter}, conventionally).

198 '''

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

200

201 polaradius = b

202

203 def _coeffs(self, auxout, auxin):

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

205 aL = self.ALorder # aka Lmax

206 if auxout == auxin:

207 return _0_0s(aL) # uncached

208

209 k = Aux._1d(auxout, auxin)

210 try: # cached

211 return AuxLat._csc[aL][k]

212 except KeyError:

213 pass

214

215 Cx = _CXcoeffs(aL)

216 try:

217 Cx = Cx[auxout][auxin]

218 except KeyError as x:

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

220

221 d = x = n = self._n

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

223 x = self._n2

224

225 def _m(m):

226 return m // 2

227 else:

228 def _m(m): # PYCHOK expected

229 return m

230

231 i = 0

232 cs = []

233 _c = cs.append

234 _p = _polynomial

235 for r in _reverange(aL):

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

237 _c(_p(x, Cx, i, j) * d)

238 d *= n

239 i = j

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

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

242 return cs

243

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

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

246

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

248

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

250 '''

251 _xinstanceof(AuxAngle, Phi=Phi)

252 # assert Phi._AUX == Aux.PHI

253 tphi = tchi = fabs(Phi.tan)

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

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

256 scsig = _sc(sig)

257 scphi = _sc(tphi)

258 if self.f > 0:

259 # The general expression for tchi is

260 # tphi * scsig - sig * scphi

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

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

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

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

265 sigtphi = sig / tphi

266 if sig < (tphi * _0_5):

267 t = tphi - sig

268 else:

269 def _asinh_2(x):

270 return asinh(x) * _0_5

271 # Still possibly dangerous cancellation in tphi - sig.

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

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

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

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

276 # Turn the crank on divided differences, substitute

277 # sphi = tphi / _sc(tphi)

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

279 e = self._e

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

281 # assert em1 != 0

282 scb = self._scbeta(tphi)

283 scphib = scphi / scb # sec(phi) / sec(beta)

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

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

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

287 try:

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

289 except ValueError: # Fsum(NAN) exception

290 Dg = _sum((atphi, atphib, t, e * t))

291 e *= atphib

292 t = atphi - e

293 if t: # sinh(0) == 0

294 Dg *= sinh(t) / t * cosh(atphi + e) * em1

295 t = float(Dg) # tphi - sig

296 tchi = _over(t * (_1_0 + sigtphi),

297 scsig + scphi * sigtphi) if t else _0_0

298 else:

299 tchi = tphi * scsig - sig * scphi

300

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

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

303

304 if _xkwds_get(diff_name, diff=False):

305 if isinf(tphi): # PYCHOK np cover

306 d = self._conformal_diff

307 else:

308 d = self.Parametric(Phi)._x_normalized

309 if d:

310 d = _over(d, Chi._x_normalized) * \

311 _over(d, Phi._x_normalized) * self._e2m1

312 Chi._diff = d

313 # assrt Chi._AUX == Aux.CHI

314 return Chi

315

316 @Property_RO

317 def _conformal_diff(self): # PYCHOK no cover

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

319 '''

320 e = self._e

321 if self.f > 0:

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

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

324 elif e:

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

326 d = _sc(ss) - ss

327 else:

328 d = _1_0

329 return d

330

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

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

333 Z = d = Zeta_d

334 if isinstance(Z, AuxAngle):

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

336 if auxin == auxout:

337 pass

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

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

340 elif exact:

341 p = Aux.power(auxout, auxin)

342 if p is None:

343 P = self._fromAux(Z) # Phi

344 Z = self._toAux(auxout, P)

345 Z._iter = P.iteration

346 else:

347 y, x = Z._yx

348 if p:

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

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

351 else:

352 cs = self._coeffs(auxout, auxin)

353 yx = Z._yx_normalized

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

355 # assert Z._yx == yx

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

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

358 # assert Z._AUX == auxout

359 return Z

360

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

362 Z = AuxPhi.fromDegrees(d)

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

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

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

366

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

368

369 def _Dq(self, tphi):

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

371 sphi = _sn(tphi)

372 if tphi > 0:

373 scphi = _sc(tphi)

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

375 if d:

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

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

378 # (1 + e2 * sphi) / ((1 - e2 * sphi * sphi) * e2m1)

379 # with atanh(e * d / (1 - e2 * sphi)) =

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

381 e2m1, ed = self._e2m1, (self._e * d)

382 if e2m1 and ed:

383 e2 = self._e2

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

385 scb = self._scbeta(tphi)

386 q = scphib = scphi / scb

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

388 q += asinh(self._e1 * d * scphib) / ed

389 else:

390 s2 = sphi * e2

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

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

393 else: # PYCHOK no cover

394 q = INF

395 else: # PYCHOK no cover

396 q = self._2_e2m12

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

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

399 return q

400

401 @Property_RO

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

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

404

405 @Property_RO

406 def _e1(self):

407 return sqrt(fabs(self._e12))

408

409 @Property_RO

410 def _e12(self):

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

412

413 @Property_RO

414 def _e12p1(self):

415 return _1_0 / self._e2m1

416

417 @Property_RO

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

419 return self.ellipsoid.e2

420

421 @Property_RO

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

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

424

425 @Property_RO

426 def _e2m1_sq2(self):

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

428

429 @Property_RO

430 def _2_e2m12(self):

431 return _2_0 / self._e2m1**2

432

433 @Property_RO

434 def _Ef_fRG_a2b2_PI_4(self):

435 E = self.ellipsoid

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

437

438 @Property_RO

439 def ellipsoid(self):

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

441 '''

442 return self._E

443

444 @Property_RO

445 def f(self):

446 '''Get the C{ellipsoid}'s flattening (C{scalar}).

447 '''

448 return self.ellipsoid.f

449

450 flattening = f

451

452 @Property_RO

453 def _fm1(self): # 1 - flattening

454 return self.ellipsoid.f1

455

456 def _fromAux(self, Zeta, **name):

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

458

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

460

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

462 '''

463 _xinstanceof(AuxAngle, Zeta=Zeta)

464 aux = Zeta._AUX

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

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

467 if f is None:

468 Phi = AuxPhi(NAN, name=n)

469 elif callable(f):

470 t = self._fm1

471 t *= f(t)

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

473 else: # assert isscalar(f)

474 y, x = Zeta._yx

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

476 # assert Phi._AUX == Aux.PHI

477 return Phi

478

479 @Property_RO

480 def _fromAuxCase(self):

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

482 '''

483 return {Aux.AUTHALIC: cbrt,

484 Aux.CONFORMAL: _MODS.utily._passarg,

485 Aux.GEOCENTRIC: self._e2m1,

486 Aux.GEOGRAPHIC: _1_0,

487 Aux.PARAMETRIC: self._fm1,

488 Aux.RECTIFYING: sqrt}

489

490 def Geocentric(self, Phi, **diff_name):

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

492

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

494

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

496 '''

497 _xinstanceof(AuxAngle, Phi=Phi)

498 # assert Phi._AUX == Aux.PHI

499 Theta = AuxTheta(Phi.y * self._e2m1, Phi.x,

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

501 if _xkwds_get(diff_name, diff=False):

502 Theta._diff = self._e2m1

503 return Theta

504

505 def Geodetic(self, Phi, **diff_name): # PYCHOK no cover

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

507

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

509

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

511 '''

512 _xinstanceof(AuxAngle, Phi=Phi)

513 # assert Phi._AUX == Aux.PHI

514 return AuxPhi(Phi, name=_xkwds_get(diff_name, name=Phi.name))

515

516 @Property_RO

517 def _n(self): # 3rd flattening

518 return self.ellipsoid.n

519

520 @Property_RO

521 def _n2(self):

522 return self._n**2

523

524 def Parametric(self, Phi, **diff_name):

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

526

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

528

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

530 '''

531 _xinstanceof(AuxAngle, Phi=Phi)

532 # assert Phi._AUX == Aux.PHI

533 Beta = AuxBeta(Phi.y * self._fm1, Phi.x,

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

535 if _xkwds_get(diff_name, diff=False):

536 Beta._diff = self._fm1

537 return Beta

538

539 Reduced = Parametric

540

541 @Property_RO

542 def _q(self): # constant _q

543 q, f = self._e12p1, self.f

544 if f:

545 e = self._e

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

547 else:

548 q += _1_0

549 return q

550

551 def _qf(self, tphi):

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

553 scb = self._scbeta(tphi)

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

555

556 def _qIntegrand(self, beta):

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

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

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

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

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

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

563 # then the integral is

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

565 # where

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

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

568 Beta = AuxBeta.fromRadians(beta)

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

570 Ax, _cv = Aux, self.convert

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

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

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

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

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

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

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

578 r = _0_0 # XXX 0 avoids NAN summation exceptions

579 return r

580

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

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

583

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

585

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

587 '''

588 Beta = self.Parametric(Phi)

589 # assert Beta._AUX == Aux.BETA

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

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

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

593 if self.f < 0:

594 a, b = b, a

595 ka, kb = kb, ka

596 ka1, kb1 = kb1, ka1

597 sb, cb = cb, sb

598 # now a, b = larger/smaller semiaxis

599 # Beta measured from larger semiaxis

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

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

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

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

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

605 sb2 = sb**2

606 cb2 = cb**2

607 da2 = ka1 + ka * sb2

608 db2 = _1_0 - kb * sb2

609 # DLMF Eq. 19.25.9

610 my = b * sb * _Ef._RFRD(cb2, db2, _1_0, kb * sb2)

611 # DLMF Eq. 19.25.10 with complementary angles

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

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

614 ka * sb / sqrt(da2))

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

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

617 if self.f < 0:

618 a, b = b, a

619 my, mx = mx, my

620 mr = (my + mx) / PI_2

621 if mr:

622 my = sin(my / mr)

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

624 else: # zero Mu

625 my, mx = _0_0, _1_0

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

627 Mu = AuxMu(my, mx, # normalized

628 name=n).copyquadrant(Phi)

629

630 if _xkwds_get(diff_name, diff=False):

631 d, x = _0_0, Beta._x_normalized

632 if x and mr:

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

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

635 * (x / Phi._x_normalized)

636 else:

637 d = mr / a

638 Mu._diff = self._fm1 * d

639 return Mu

640

641 def RectifyingRadius(self, exact=False):

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

643

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

645 otherwise the I{Taylor} series.

646

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

648 as the ellipsoid axes).

649 '''

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

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

652

653 @Property_RO

654 def _RectifyingR(self):

655 m = self.ALorder

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

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

658

659 def _scbeta(self, tphi):

660 return _sc(self._fm1 * tphi)

661

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

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

664

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

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

667

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

669 '''

670 _xinstanceof(AuxAngle, Phi=Phi)

671 # assert Phi._AUX == Aux.PHI

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

673 m = _toAuxCase.get(auxout, None)

674 if m: # callable

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

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

677 A = AuxPhi(Phi, name=n)

678 else: # auxout?

679 A = AuxPhi(NAN, name=n)

680 # assert A._AUX == auxout

681 return A

682

683 def _toZeta(self, zetaux):

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

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

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

687 '''

688 class _AuxPhy(AuxPhi):

689 # lean C{AuxPhi} instance.

690 # _diff = _1_0

691 # _x = _1_0

692

693 def __init__(self, tphi): # PYCHOK signature

694 self._y = tphi

695

696 m = _toAuxCase.get(zetaux, None)

697 if m: # callable

698

699 def _toZeta(tphi):

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

701

702 elif zetaux == Aux.GEODETIC: # GEOGRAPHIC

703 _toZeta = _AuxPhy

704

705 else: # zetaux?

706

707 def _toZeta(unused): # PYCHOK expected

708 return _AuxPhy(NAN)

709

710 return _toZeta

711

712

713# switch(auxout): ...

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

715 Aux.CONFORMAL: AuxLat.Conformal,

716 Aux.GEOCENTRIC: AuxLat.Geocentric,

717 Aux.PARAMETRIC: AuxLat.Parametric,

718 Aux.RECTIFYING: AuxLat.Rectifying}

719

720

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

722 sz, cz = Zeta._yx # isnormal

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

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

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

726 if isfinite(x):

727 U0, U1 = Fsum(), Fsum()

728 # assert len(cs) == K

729 for r in _reverange(K):

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

731 U1, U0 = U0, -U1

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

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

734 # fm1 = 0 if sinp else 1

735 if sinp:

736 U0 *= sz * cz * _2_0

737 else:

738 U0 *= x * _0_5

739 U0 -= U1

740 x = float(U0)

741 return x

742

743

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

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

746 '''

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

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

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

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

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

752 # assert _coeffs.ALorder == aL

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

754 return _coeffs

755

756

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

758 # Newton's method fro AuxLat._fromAux

759 try:

760 _lg2 = _log2

761 _abs = fabs

762 tz = _abs(Zeta.tan)

763 tphi = tz / tphi # **)

764 ltz = _lg2(tz) # **)

765 ltphi = _lg2(tphi) # **)

766 ltmin = min(ltphi, MIN_EXP)

767 ltmax = max(ltphi, MAX_EXP)

768# auxin = Zeta._AUX

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

770 _TOL, _xp2 = _EPSqrt, _exp2

771 for i in range(1, _TRIPS): # up to 1 Ki!

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

773 Z = _toZeta(tphi)

774 # assert Z._AUX == auxin

775 t, s_ = Z.tan, s

776 if t > tz:

777 ltmax, s = ltphi, +1

778 elif t < tz:

779 ltmin, s = ltphi, -1

780 else:

781 break

782 # derivative dtan(Z)/dtan(Phi)

783 # to dlog(tan(Z))/dlog(tan(Phi))

784 d = (ltz - _lg2(t)) * t # **)

785 if d:

786 d = d / (Z.diff * tphi) # **)

787 ltphi += d

788 tphi = _xp2(ltphi)

789 if _abs(d) < _TOL:

790 i += 1

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

792 Z = _toZeta(tphi)

793 tphi -= _over(Z.tan - tz, Z.diff)

794 break

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

796 s, n = 0, (i + 2)

797 ltphi = (ltmin + ltmax) * __2

798 tphi = _xp2(ltphi)

799 else:

800 i = _TRIPS

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

802 Phi._iter = i

803 except (ValueError, ZeroDivisionError): # **) zero t, tphi, tz or Z.diff

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

805 Phi._iter = 0

806 # assert Phi._AUX == Aux.PHI

807 return Phi

808

809

810_f, _u = float, _Ufloats()

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

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

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

814 4 / _f(15), _1__f3),

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

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

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

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

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

820 del _f, _u, _Ufloats, _1__f3

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

822

823# AuxLat._Lmax = max(_AR2Coeffs.keys()) # == max(ALorder)

824

825__all__ += _ALL_DOCS(AuxLat)

826

827# **) MIT License

828#

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

830#

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

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

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

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

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

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

837#

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

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

840#

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

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

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

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

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

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

847# OTHER DEALINGS IN THE SOFTWARE.