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, atan1

21 from pygeodesy.basics import _passarg, _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, _Fsumf_, _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 _isDegrees, _isRadius, Degrees, Meter

36# from pygeodesy.utily import atan1 # from .auxily

37

38 from math import asinh, 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__ = '24.04.14'

48

49_TRIPS = 1024 # XXX 2 or 3?

50

51

52 class AuxLat(AuxAngle):

53 '''Base class for accurate conversion between I{Auxiliary} latitudes

54 on an ellipsoid.

55

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

57 maintain precision near the poles, I{Authalic} latitude C{Xi},

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

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

60

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

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

63 '''

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

65 _E = _WGS84.ellipsoid

66# _Lmax = 0 # overwritten below

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 _isRadius(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 C{ellipsoid}'s 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 = _over(atan1(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 = _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 C{ellipsoid}'s 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 if auxout == auxin:

208 return _0_0s(aL) # uncached

209

210 k = Aux._1d(auxout, auxin)

211 try: # cached

212 return AuxLat._csc[aL][k]

213 except KeyError:

214 pass

215

216 Cx = _CXcoeffs(aL)

217 try:

218 Cx = Cx[auxout][auxin]

219 except KeyError as x:

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

221

222 d = x = n = self._n

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

224 x = self._n2

225

226 def _m(aL):

227 for m in _reverange(aL):

228 yield m // 2

229 else:

230 _m = _reverange # PYCHOK expected

231

232 i = 0

233 cs = []

234 _c = cs.append

235 _p = _polynomial

236 for m in _m(aL):

237 j = i + m + 1 # order m = j - i - 1

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

239 d *= n

240 i = j

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

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

243 return cs

244

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

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

247

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

249

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

251 '''

252 _xinstanceof(AuxAngle, Phi=Phi)

253 # assert Phi._AUX == Aux.PHI

254 tphi = tchi = fabs(Phi.tan)

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

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

257 scsig = _sc(sig)

258 scphi = _sc(tphi)

259 if self.f > 0:

260 # The general expression for tchi is

261 # tphi * scsig - sig * scphi

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

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

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

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

266 sigtphi = sig / tphi

267 if sig < (tphi * _0_5):

268 t = tphi - sig

269 else:

270 def _asinh_2(x):

271 return asinh(x) * _0_5

272 # Still possibly dangerous cancellation in tphi - sig.

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

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

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

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

277 # Turn the crank on divided differences, substitute

278 # sphi = tphi / _sc(tphi)

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

280 e = self._e

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

282 # assert em1 != 0

283 scb = self._scbeta(tphi)

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

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

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

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

288 try:

289 Dg = _Fsumf_(atphi, atphib, t, e * t)

290 except ValueError: # Fsum(NAN) exception

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

292 e *= atphib

293 t = atphi - e

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

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

296 t = float(Dg) # tphi - sig

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

298 scsig + scphi * sigtphi) if t else _0_0

299 else:

300 tchi = tphi * scsig - sig * scphi

301

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

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

304

305 if _xkwds_get(diff_name, diff=False):

306 if isinf(tphi): # PYCHOK np cover

307 d = self._conformal_diff

308 else:

309 d = self.Parametric(Phi)._x_normalized

310 if d:

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

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

313 Chi._diff = d

314 # assrt Chi._AUX == Aux.CHI

315 return Chi

316

317 @Property_RO

318 def _conformal_diff(self): # PYCHOK no cover

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

320 '''

321 e = self._e

322 if self.f > 0:

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

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

325 elif e:

326 ss = sinh(-atan1(e) * e)

327 d = _sc(ss) - ss

328 else:

329 d = _1_0

330 return d

331

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

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

334 Z = d = Zeta_d

335 if isinstance(Z, AuxAngle):

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

337 if auxin == auxout:

338 pass

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

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

341 elif exact:

342 p = Aux.power(auxout, auxin)

343 if p is None:

344 P = self._fromAux(Z) # Phi

345 Z = self._toAux(auxout, P)

346 Z._iter = P.iteration

347 else:

348 y, x = Z._yx

349 if p:

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

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

352 else:

353 cs = self._coeffs(auxout, auxin)

354 yx = Z._yx_normalized

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

356 # assert Z._yx == yx

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

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

359 # assert Z._AUX == auxout

360 return Z

361

362 elif _isDegrees(d):

363 Z = AuxPhi.fromDegrees(d)

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

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

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

367

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

369

370 def _Dq(self, tphi):

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

372 sphi = _sn(tphi)

373 if tphi > 0:

374 scphi = _sc(tphi)

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

376 if d:

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

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

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

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

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

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

383 if e2m1 and ed:

384 e2 = self._e2

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

386 scb = self._scbeta(tphi)

387 q = scphib = scphi / scb

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

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

390 else:

391 s2 = sphi * e2

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

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

394 else: # PYCHOK no cover

395 q = INF

396 else: # PYCHOK no cover

397 q = self._2_e2m12

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

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

400 return q

401

402 @Property_RO

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

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

405

406 @Property_RO

407 def _e1(self):

408 return sqrt(fabs(self._e12))

409

410 @Property_RO

411 def _e12(self):

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

413

414 @Property_RO

415 def _e12p1(self):

416 return _1_0 / self._e2m1

417

418 @Property_RO

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

420 return self.ellipsoid.e2

421

422 @Property_RO

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

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

425

426 @Property_RO

427 def _e2m1_sq2(self):

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

429

430 @Property_RO

431 def _2_e2m12(self):

432 return _2_0 / self._e2m1**2

433

434 @Property_RO

435 def _Ef_fRG_a2b2_PI_4(self):

436 E = self.ellipsoid

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

438

439 @Property_RO

440 def ellipsoid(self):

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

442 '''

443 return self._E

444

445 @Property_RO

446 def f(self):

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

448 '''

449 return self.ellipsoid.f

450

451 flattening = f

452

453 @Property_RO

454 def _fm1(self): # 1 - flattening

455 return self.ellipsoid.f1

456

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

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

459

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

461

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

463 '''

464 _xinstanceof(AuxAngle, Zeta=Zeta)

465 aux = Zeta._AUX

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

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

468 if f is None:

469 Phi = AuxPhi(NAN, name=n)

470 elif callable(f):

471 t = self._fm1

472 t *= f(t)

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

474 else: # assert isscalar(f)

475 y, x = Zeta._yx

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

477 # assert Phi._AUX == Aux.PHI

478 return Phi

479

480 @Property_RO

481 def _fromAuxCase(self):

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

483 '''

484 return {Aux.AUTHALIC: cbrt,

485 Aux.CONFORMAL: _passarg,

486 Aux.GEOCENTRIC: self._e2m1,

487 Aux.GEOGRAPHIC: _1_0,

488 Aux.PARAMETRIC: self._fm1,

489 Aux.RECTIFYING: sqrt}

490

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

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

493

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

495

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

497 '''

498 _xinstanceof(AuxAngle, Phi=Phi)

499 # assert Phi._AUX == Aux.PHI

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

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

502 if _xkwds_get(diff_name, diff=False):

503 Theta._diff = self._e2m1

504 return Theta

505

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

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

508

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

510

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

512 '''

513 _xinstanceof(AuxAngle, Phi=Phi)

514 # assert Phi._AUX == Aux.PHI

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

516

517 @Property_RO

518 def _n(self): # 3rd flattening

519 return self.ellipsoid.n

520

521 @Property_RO

522 def _n2(self):

523 return self._n**2

524

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

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

527

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

529

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

531 '''

532 _xinstanceof(AuxAngle, Phi=Phi)

533 # assert Phi._AUX == Aux.PHI

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

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

536 if _xkwds_get(diff_name, diff=False):

537 Beta._diff = self._fm1

538 return Beta

539

540 Reduced = Parametric

541

542 @Property_RO

543 def _q(self): # constant _q

544 q, f = self._e12p1, self.f

545 if f:

546 e = self._e

547 q += _over(asinh(self._e1) if f > 0 else atan1(e), e)

548 else:

549 q += _1_0

550 return q

551

552 def _qf(self, tphi):

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

554 scb = self._scbeta(tphi)

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

556

557 def _qIntegrand(self, beta):

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

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

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

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

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

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

564 # then the integral is

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

566 # where

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

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

569 Beta = AuxBeta.fromRadians(beta)

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

571 Ax, _cv = Aux, self.convert

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

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

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

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

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

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

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

579 r = _0_0 # XXX 0 avoids NAN summation exceptions

580 return r

581

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

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

584

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

586

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

588 '''

589 Beta = self.Parametric(Phi)

590 # assert Beta._AUX == Aux.BETA

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

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

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

594 if self.f < 0:

595 a, b = b, a

596 ka, kb = kb, ka

597 ka1, kb1 = kb1, ka1

598 sb, cb = cb, sb

599 # now a, b = larger/smaller semiaxis

600 # Beta measured from larger semiaxis

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

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

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

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

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

606 sb2 = sb**2

607 cb2 = cb**2

608 da2 = ka1 + ka * sb2

609 db2 = _1_0 - kb * sb2

610 # DLMF Eq. 19.25.9

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

612 # DLMF Eq. 19.25.10 with complementary angles

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

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

615 ka * sb / sqrt(da2))

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

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

618 if self.f < 0:

619 a, b = b, a

620 my, mx = mx, my

621 mr = (my + mx) / PI_2

622 if mr:

623 my = sin(my / mr)

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

625 else: # zero Mu

626 my, mx = _0_0, _1_0

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

628 Mu = AuxMu(my, mx, # normalized

629 name=n).copyquadrant(Phi)

630

631 if _xkwds_get(diff_name, diff=False):

632 d, x = _0_0, Beta._x_normalized

633 if x and mr:

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

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

636 * (x / Phi._x_normalized)

637 else:

638 d = mr / a

639 Mu._diff = self._fm1 * d

640 return Mu

641

642 def RectifyingRadius(self, exact=False):

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

644

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

646 otherwise the I{Taylor} series.

647

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

649 as the ellipsoid axes).

650 '''

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

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

653

654 @Property_RO

655 def _RectifyingR(self):

656 m = self.ALorder

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

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

659

660 def _scbeta(self, tphi):

661 return _sc(self._fm1 * tphi)

662

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

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

665

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

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

668

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

670 '''

671 _xinstanceof(AuxAngle, Phi=Phi)

672 # assert Phi._AUX == Aux.PHI

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

674 m = _toAuxCase.get(auxout, None)

675 if m: # callable

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

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

678 A = AuxPhi(Phi, name=n)

679 else: # auxout?

680 A = AuxPhi(NAN, name=n)

681 # assert A._AUX == auxout

682 return A

683

684 def _toZeta(self, zetaux):

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

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

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

688 '''

689 class _AuxPhy(AuxPhi):

690 # lean C{AuxPhi} instance.

691 # _diff = _1_0

692 # _x = _1_0

693

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

695 self._y = tphi

696

697 m = _toAuxCase.get(zetaux, None)

698 if m: # callable

699

700 def _toZeta(tphi):

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

702

703 elif zetaux == Aux.GEODETIC: # GEOGRAPHIC

704 _toZeta = _AuxPhy

705

706 else: # zetaux?

707

708 def _toZeta(unused): # PYCHOK expected

709 return _AuxPhy(NAN)

710

711 return _toZeta

712

713

714# switch(auxout): ...

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

716 Aux.CONFORMAL: AuxLat.Conformal,

717 Aux.GEOCENTRIC: AuxLat.Geocentric,

718 Aux.PARAMETRIC: AuxLat.Parametric,

719 Aux.RECTIFYING: AuxLat.Rectifying}

720

721

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

723 sz, cz = Zeta._yx # isnormal

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

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

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

727 if isfinite(x):

728 U0, U1 = Fsum(), Fsum()

729 # assert len(cs) == K

730 for r in _reverange(K):

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

732 U1, U0 = U0, -U1

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

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

735 # fm1 = 0 if sinp else 1

736 if sinp:

737 U0 *= sz * cz * _2_0

738 else:

739 U0 *= x * _0_5

740 U0 -= U1

741 x = float(U0)

742 return x

743

744

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

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

747 '''

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

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

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

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

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

753 # assert _coeffs.ALorder == aL

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

755 return _coeffs

756

757

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

759 # Newton's method fro AuxLat._fromAux

760 try:

761 _lg2 = _log2

762 _abs = fabs

763 tz = _abs(Zeta.tan)

764 tphi = tz / tphi # **)

765 ltz = _lg2(tz) # **)

766 ltphi = _lg2(tphi) # **)

767 ltmin = min(ltphi, MIN_EXP)

768 ltmax = max(ltphi, MAX_EXP)

769# auxin = Zeta._AUX

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

771 _TOL, _xp2 = _EPSqrt, _exp2

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

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

774 Z = _toZeta(tphi)

775 # assert Z._AUX == auxin

776 t, s_ = Z.tan, s

777 if t > tz:

778 ltmax, s = ltphi, +1

779 elif t < tz:

780 ltmin, s = ltphi, -1

781 else:

782 break

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

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

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

786 if d:

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

788 ltphi += d

789 tphi = _xp2(ltphi)

790 if _abs(d) < _TOL:

791 i += 1

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

793 Z = _toZeta(tphi)

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

795 break

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

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

798 ltphi = (ltmin + ltmax) * __2

799 tphi = _xp2(ltphi)

800 else:

801 i = _TRIPS

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

803 Phi._iter = i

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

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

806 Phi._iter = 0

807 # assert Phi._AUX == Aux.PHI

808 return Phi

809

810

811_f, _u = float, _Ufloats()

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

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

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

815 4 / _f(15), _1__f3),

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

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

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

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

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

821 del _f, _u, _Ufloats, _1__f3

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

823

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

825

826__all__ += _ALL_DOCS(AuxLat)

827

828# **) MIT License

829#

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

831#

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

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

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

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

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

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

838#

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

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

841#

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

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

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

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

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

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

848# OTHER DEALINGS IN THE SOFTWARE.