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 _reverange, _xinstanceof, _passarg

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

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

27 from pygeodesy.elliptic import Elliptic as _Ef

28 from pygeodesy.errors import AuxError, _xkwds_not, _xkwds_pop2, _Xorder

29# from pygeodesy.fmath import cbrt # from .karney

30 from pygeodesy.fsums import Fsum, _Fsumf_, _sum

31# from pygeodesy.internals import _passarg # from .basics

32 from pygeodesy.interns import NN, _DOT_, _not_scalar_, _UNDER_

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

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

35# from pygeodesy.named import _name__ # from .datums

36 from pygeodesy.props import Property, Property_RO, _update_all

37 from pygeodesy.units import _isDegrees, _isRadius, Degrees, Meter

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

39

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

41 try:

42 from math import exp2 as _exp2

43 except ImportError: # Python 3.11-

44

45 def _exp2(x):

46 return pow(_2_0, x)

47

48__all__ = ()

49__version__ = '24.05.24'

50

51_TRIPS = 1024 # XXX 2 or 3?

52

53

54 class AuxLat(AuxAngle):

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

56 on an ellipsoid.

57

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

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

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

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

62

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

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

65 '''

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

67 _E = _WGS84.ellipsoid

68# _Lmax = 0 # overwritten below

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

70

71 def __init__(self, a_earth=_WGS84, f=None, b=None, **ALorder_name):

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

73

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

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

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

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

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

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

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

81 is not scalar.

82 @kwarg ALorder_name: Optional C{B{name}=NN} (C{str}) and optional

83 keyword argument C{B{ALorder}=6} for the order of

84 this L{AuxLat}, see property C{ALorder}.

85 '''

86 if ALorder_name:

87 M = self._mAL

88 m, name = _xkwds_pop2(ALorder_name, ALorder=M)

89 if m != M:

90 self.ALorder = m

91 else:

92 name = NN

93 try:

94 if a_earth is not _WGS84:

95 n = _name__(name, name__=AuxLat)

96 if b is f is None:

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

98 elif _isRadius(a_earth):

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

100 else:

101 raise ValueError(_not_scalar_)

102 self._E = E

103 elif not (b is f is None):

104 # _UnexpectedError into AuxError

105 name = _name__(name, **_xkwds_not(None, b=b, f=f))

106

107 if name:

108 self.name = name

109 except Exception as x:

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

111

112 @Property_RO

113 def a(self):

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

115 '''

116 return self.ellipsoid.a

117

118 equatoradius = a

119

120 @Property

121 def ALorder(self):

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

123 '''

124 return self._mAL

125

126 @ALorder.setter # PYCHOK setter!

127 def ALorder(self, order):

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

129 '''

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

131 if self._mAL != m:

132 _update_all(self)

133 self._mAL = m

134

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

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

137 f = self.f

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

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

140 e = self._e

141 s = _over(atan1(e * s) if f < 0 else asinh(self._e1 * s), e)

142 return s

143

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

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

146

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

148 @kwarg diff_name: Use C{B{diff}=True} to set C{diff}

149 optional C{B{name}=NN}.

150

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

152 '''

153 _xinstanceof(AuxAngle, Phi=Phi)

154 # assert Phi._AUX == Aux.PHI

155 tphi = fabs(Phi.tan)

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

157 y, x = Phi._yx_normalized

158 q = self._q

159 qv = self._qf(tphi)

160 Dq2 = self._Dq(tphi)

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

162 d, n = _diff_name2(Phi, **diff_name)

163 Xi = AuxXi(copysign(qv, Phi.y), x * sqrt(Dq2), name=n)

164

165 if d:

166 if isnan(tphi):

167 d = self._e2m1_sq2

168 else:

169 c = self.Parametric(Phi)._x_normalized

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

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

172 Xi._diff = d

173 else:

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

175 # assert Xi._AUX == Aux.XI

176 return Xi

177

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

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

180

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

182 the I{Taylor} series.

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

184

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

186 units as the ellipsoid axes).

187

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

189 exceeds C{f_max}.

190 '''

191 f = self.f

192 if exact or not f:

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

194 elif fabs(f) < f_max:

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

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

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

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

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

200 m = self.ALorder

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

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

203 else:

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

205 return c2 * _0_5

206

207 @Property_RO

208 def b(self):

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

210 '''

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

212

213 polaradius = b

214

215 def _coeffs(self, auxout, auxin):

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

217 aL = self.ALorder # aka Lmax

218 if auxout == auxin:

219 return _0_0s(aL) # uncached

220

221 k = Aux._1d(auxout, auxin)

222 try: # cached

223 return AuxLat._csc[aL][k]

224 except KeyError:

225 pass

226

227 Cx = _CXcoeffs(aL)

228 try:

229 Cx = Cx[auxout][auxin]

230 except KeyError as x:

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

232

233 d = x = n = self._n

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

235 x = self._n2

236

237 def _m(aL):

238 for m in _reverange(aL):

239 yield m // 2

240 else:

241 _m = _reverange # PYCHOK expected

242

243 i = 0

244 cs = []

245 _c = cs.append

246 _p = _polynomial

247 for m in _m(aL):

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

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

250 d *= n

251 i = j

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

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

254 return cs

255

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

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

258

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

260 @kwarg diff_name: Use C{B{diff}=True} to set C{diff}

261 and an optional C{B{name}=NN}.

262

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

264 '''

265 _xinstanceof(AuxAngle, Phi=Phi)

266 # assert Phi._AUX == Aux.PHI

267 tphi = tchi = fabs(Phi.tan)

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

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

270 scsig = _sc(sig)

271 scphi = _sc(tphi)

272 if self.f > 0:

273 # The general expression for tchi is

274 # tphi * scsig - sig * scphi

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

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

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

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

279 sigtphi = sig / tphi

280 if sig < (tphi * _0_5):

281 t = tphi - sig

282 else:

283 def _asinh_2(x):

284 return asinh(x) * _0_5

285 # Still possibly dangerous cancellation in tphi - sig.

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

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

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

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

290 # Turn the crank on divided differences, substitute

291 # sphi = tphi / _sc(tphi)

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

293 e = self._e

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

295 # assert em1 != 0

296 scb = self._scbeta(tphi)

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

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

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

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

301 try:

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

303 except ValueError: # Fsum(NAN) exception

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

305 e *= atphib

306 t = atphi - e

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

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

309 t = float(Dg) # tphi - sig

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

311 scsig + scphi * sigtphi) if t else _0_0

312 else:

313 tchi = tphi * scsig - sig * scphi

314

315 d, n = _diff_name2(Phi, **diff_name)

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

317

318 if d:

319 if isinf(tphi): # PYCHOK np cover

320 d = self._conformal_diff

321 else:

322 d = self.Parametric(Phi)._x_normalized

323 if d:

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

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

326 Chi._diff = d

327 # assrt Chi._AUX == Aux.CHI

328 return Chi

329

330 @Property_RO

331 def _conformal_diff(self): # PYCHOK no cover

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

333 '''

334 e = self._e

335 if self.f > 0:

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

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

338 elif e:

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

340 d = _sc(ss) - ss

341 else:

342 d = _1_0

343 return d

344

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

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

347 Z = d = Zeta_d

348 if isinstance(Z, AuxAngle):

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

350 if auxin == auxout:

351 pass

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

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

354 elif exact:

355 p = Aux.power(auxout, auxin)

356 if p is None:

357 P = self._fromAux(Z) # Phi

358 Z = self._toAux(auxout, P)

359 Z._iter = P.iteration

360 else:

361 y, x = Z._yx

362 if p:

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

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

365 else:

366 cs = self._coeffs(auxout, auxin)

367 yx = Z._yx_normalized

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

369 # assert Z._yx == yx

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

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

372 # assert Z._AUX == auxout

373 return Z

374

375 elif _isDegrees(d):

376 Z = AuxPhi.fromDegrees(d)

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

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

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

380

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

382

383 def _Dq(self, tphi):

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

385 sphi = _sn(tphi)

386 if tphi > 0:

387 scphi = _sc(tphi)

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

389 if d:

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

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

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

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

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

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

396 if e2m1 and ed:

397 e2 = self._e2

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

399 scb = self._scbeta(tphi)

400 q = scphib = scphi / scb

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

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

403 else:

404 s2 = sphi * e2

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

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

407 else: # PYCHOK no cover

408 q = INF

409 else: # PYCHOK no cover

410 q = self._2_e2m12

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

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

413 return q

414

415 @Property_RO

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

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

418

419 @Property_RO

420 def _e1(self):

421 return sqrt(fabs(self._e12))

422

423 @Property_RO

424 def _e12(self):

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

426

427 @Property_RO

428 def _e12p1(self):

429 return _1_0 / self._e2m1

430

431 @Property_RO

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

433 return self.ellipsoid.e2

434

435 @Property_RO

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

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

438

439 @Property_RO

440 def _e2m1_sq2(self):

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

442

443 @Property_RO

444 def _2_e2m12(self):

445 return _2_0 / self._e2m1**2

446

447 @Property_RO

448 def _Ef_fRG_a2b2_PI_4(self):

449 E = self.ellipsoid

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

451

452 @Property_RO

453 def ellipsoid(self):

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

455 '''

456 return self._E

457

458 @Property_RO

459 def f(self):

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

461 '''

462 return self.ellipsoid.f

463

464 flattening = f

465

466 @Property_RO

467 def _fm1(self): # 1 - flattening

468 return self.ellipsoid.f1

469

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

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

472

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

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

475

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

477 '''

478 _xinstanceof(AuxAngle, Zeta=Zeta)

479 aux = Zeta._AUX

480 n = _name__(name, _or_nameof=Zeta)

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

482 if f is None:

483 Phi = AuxPhi(NAN, name=n)

484 elif callable(f):

485 t = self._fm1

486 t *= f(t)

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

488 else: # assert isscalar(f)

489 y, x = Zeta._yx

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

491 # assert Phi._AUX == Aux.PHI

492 return Phi

493

494 @Property_RO

495 def _fromAuxCase(self):

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

497 '''

498 return {Aux.AUTHALIC: cbrt,

499 Aux.CONFORMAL: _passarg,

500 Aux.GEOCENTRIC: self._e2m1,

501 Aux.GEOGRAPHIC: _1_0,

502 Aux.PARAMETRIC: self._fm1,

503 Aux.RECTIFYING: sqrt}

504

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

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

507

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

509 @kwarg diff_name: Use C{B{diff}=True} to set C{diff}

510 and an optional C{B{name}=NN}.

511

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

513 '''

514 _xinstanceof(AuxAngle, Phi=Phi)

515 # assert Phi._AUX == Aux.PHI

516 d, n = _diff_name2(Phi, **diff_name)

517 Theta = AuxTheta(Phi.y * self._e2m1, Phi.x, name=n)

518 if d:

519 Theta._diff = self._e2m1

520 return Theta

521

522 def Geodetic(self, Phi, **name): # PYCHOK no cover

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

524

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

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

527

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

529 '''

530 _xinstanceof(AuxAngle, Phi=Phi)

531 # assert Phi._AUX == Aux.PHI

532 _, n = _diff_name2(Phi, **name)

533 return AuxPhi(Phi, name=n)

534

535 @Property_RO

536 def _n(self): # 3rd flattening

537 return self.ellipsoid.n

538

539 @Property_RO

540 def _n2(self):

541 return self._n**2

542

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

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

545

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

547 @kwarg diff_name: Use C{B{diff}=True} to set C{diff}

548 and an optional C{B{name}=NN}.

549

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

551 '''

552 _xinstanceof(AuxAngle, Phi=Phi)

553 # assert Phi._AUX == Aux.PHI

554 d, n = _diff_name2(Phi, **diff_name)

555 Beta = AuxBeta(Phi.y * self._fm1, Phi.x, name=n)

556 if d:

557 Beta._diff = self._fm1

558 return Beta

559

560 Reduced = Parametric

561

562 @Property_RO

563 def _q(self): # constant _q

564 q, f = self._e12p1, self.f

565 if f:

566 e = self._e

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

568 else:

569 q += _1_0

570 return q

571

572 def _qf(self, tphi):

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

574 scb = self._scbeta(tphi)

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

576

577 def _qIntegrand(self, beta):

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

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

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

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

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

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

584 # then the integral is

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

586 # where

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

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

589 Beta = AuxBeta.fromRadians(beta)

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

591 Ax, _cv = Aux, self.convert

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

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

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

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

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

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

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

599 r = _0_0 # XXX 0 avoids NAN summation exceptions

600 return r

601

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

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

604

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

606 @kwarg diff_name: Use C{B{diff}=True} to set C{diff}

607 and an optional C{B{name}=NN}.

608

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

610 '''

611 Beta = self.Parametric(Phi)

612 # assert Beta._AUX == Aux.BETA

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

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

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

616 if self.f < 0:

617 a, b = b, a

618 ka, kb = kb, ka

619 ka1, kb1 = kb1, ka1

620 sb, cb = cb, sb

621 # now a, b = larger/smaller semiaxis

622 # Beta measured from larger semiaxis

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

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

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

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

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

628 sb2 = sb**2

629 cb2 = cb**2

630 da2 = ka1 + ka * sb2

631 db2 = _1_0 - kb * sb2

632 # DLMF Eq. 19.25.9

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

634 # DLMF Eq. 19.25.10 with complementary angles

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

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

637 ka * sb / sqrt(da2))

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

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

640 if self.f < 0:

641 a, b = b, a

642 my, mx = mx, my

643 mr = (my + mx) / PI_2

644 if mr:

645 my = sin(my / mr)

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

647 else: # zero Mu

648 my, mx = _0_0, _1_0

649 d, n = _diff_name2(Phi, **diff_name)

650 Mu = AuxMu(my, mx, # normalized

651 name=n).copyquadrant(Phi)

652 if d:

653 d, x = _0_0, Beta._x_normalized

654 if x and mr:

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

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

657 * (x / Phi._x_normalized)

658 else:

659 d = mr / a

660 Mu._diff = self._fm1 * d

661 return Mu

662

663 def RectifyingRadius(self, exact=False):

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

665

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

667 otherwise the I{Taylor} series.

668

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

670 as the ellipsoid axes).

671 '''

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

673 return Meter(r, name__=self.RectifyingRadius)

674

675 @Property_RO

676 def _RectifyingR(self):

677 m = self.ALorder

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

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

680

681 def _scbeta(self, tphi):

682 return _sc(self._fm1 * tphi)

683

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

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

686

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

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

689 @kwarg diff_name: Use C{B{diff}=True} to set C{diff}

690 and an optional C{B{name}=NN}.

691

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

693 '''

694 _xinstanceof(AuxAngle, Phi=Phi)

695 # assert Phi._AUX == Aux.PHI

696 d, n = _diff_name2(Phi, **diff_name)

697 m = _toAuxCase.get(auxout, None)

698 if m: # callable

699 A = m(self, Phi, diff=d, name=n)

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

701 A = AuxPhi(Phi, name=n)

702 else: # auxout?

703 A = AuxPhi(NAN, name=n)

704 # assert A._AUX == auxout

705 return A

706

707 def _toZeta(self, zetaux):

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

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

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

711 '''

712 class _AuxPhy(AuxPhi):

713 # lean C{AuxPhi} instance.

714 # _diff = _1_0

715 # _x = _1_0

716

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

718 self._y = tphi

719

720 m = _toAuxCase.get(zetaux, None)

721 if m: # callable

722

723 def _toZeta(tphi):

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

725

726 elif zetaux == Aux.GEODETIC: # GEOGRAPHIC

727 _toZeta = _AuxPhy

728

729 else: # zetaux?

730

731 def _toZeta(unused): # PYCHOK expected

732 return _AuxPhy(NAN)

733

734 return _toZeta

735

736

737# switch(auxout): ...

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

739 Aux.CONFORMAL: AuxLat.Conformal,

740 Aux.GEOCENTRIC: AuxLat.Geocentric,

741 Aux.PARAMETRIC: AuxLat.Parametric,

742 Aux.RECTIFYING: AuxLat.Rectifying}

743

744

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

746 sz, cz = Zeta._yx # isnormal

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

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

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

750 if isfinite(x):

751 U0, U1 = Fsum(), Fsum()

752 # assert len(cs) == K

753 for r in _reverange(K):

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

755 U1, U0 = U0, -U1

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

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

758 # fm1 = 0 if sinp else 1

759 if sinp:

760 U0 *= sz * cz * _2_0

761 else:

762 U0 *= x * _0_5

763 U0 -= U1

764 x = float(U0)

765 return x

766

767

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

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

770 '''

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

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

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

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

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

776 # assert _coeffs.ALorder == aL

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

778 return _coeffs

779

780

781 def _diff_name2(Phi, diff=False, **name):

782 '''(INTERNAL) Get C{{Bdiff}=False} and C{B{name}=NN}.

783 '''

784 n = _name__(name, _or_nameof=Phi) # if name else Phi.name

785 return diff, n

786

787

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

789 # Newton's method from AuxLat._fromAux

790 try:

791 _lg2 = _log2

792 _abs = fabs

793 tz = _abs(Zeta.tan)

794 tphi = tz / tphi # **)

795 ltz = _lg2(tz) # **)

796 ltphi = _lg2(tphi) # **)

797 ltmin = min(ltphi, MIN_EXP)

798 ltmax = max(ltphi, MAX_EXP)

799# auxin = Zeta._AUX

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

801 _TOL, _xp2 = _EPSqrt, _exp2

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

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

804 Z = _toZeta(tphi)

805 # assert Z._AUX == auxin

806 t, s_ = Z.tan, s

807 if t > tz:

808 ltmax, s = ltphi, +1

809 elif t < tz:

810 ltmin, s = ltphi, -1

811 else:

812 break

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

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

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

816 if d:

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

818 ltphi += d

819 tphi = _xp2(ltphi)

820 if _abs(d) < _TOL:

821 i += 1

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

823 Z = _toZeta(tphi)

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

825 break

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

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

828 ltphi = (ltmin + ltmax) * __2

829 tphi = _xp2(ltphi)

830 else:

831 i = _TRIPS

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

833 Phi._iter = i

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

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

836 Phi._iter = 0

837 # assert Phi._AUX == Aux.PHI

838 return Phi

839

840

841_f, _u = float, _Ufloats()

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

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

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

845 4 / _f(15), _1__f3),

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

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

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

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

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

851 del _f, _u, _Ufloats, _1__f3

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

853

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

855

856__all__ += _ALL_DOCS(AuxLat)

857

858# **) MIT License

859#

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

861#

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

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

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

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

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

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

868#

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

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

871#

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

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

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

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

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

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

878# OTHER DEALINGS IN THE SOFTWARE.