Coverage for pygeodesy/auxilats/auxLat.py: 93%
410 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-15 09:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-15 09:43 -0400
2# -*- coding: utf-8 -*-
4u'''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>}
6in I{GeographicLib version 2.2+}.
8Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2022-2023) and licensed
9under the MIT/X11 License. For more information, see the U{GeographicLib
10<https://GeographicLib.SourceForge.io>} documentation.
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
16from __future__ import division as _; del _ # PYCHOK semicolon
18from pygeodesy.auxilats.auxAngle import AuxAngle, AuxBeta, AuxChi, _AuxClass, \
19 AuxMu, AuxPhi, AuxTheta, AuxXi
20from pygeodesy.auxilats.auxily import Aux, _sc, _sn, _Ufloats
21from pygeodesy.basics import isscalar, _reverange, _xinstanceof
22from 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
25from pygeodesy.datums import _ellipsoidal_datum, _WGS84, Ellipsoid
26# from pygeodesy.ellipsoids import Ellipsoid # from .datums
27from pygeodesy.elliptic import Elliptic as _Ef
28from pygeodesy.errors import AuxError, _xkwds, _xkwds_get, _Xorder
29# from pygeodesy.fmath import cbrt # from .karney
30from pygeodesy.fsums import Fsum, _sum
31from pygeodesy.karney import _2cos2x, _ALL_DOCS, cbrt, _MODS
32from pygeodesy.interns import NN, _DOT_, _UNDER_ # _earth_
33# from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS # from .karney
34from pygeodesy.props import Property, Property_RO, _update_all
35from pygeodesy.units import Degrees, Meter
36# from pygeodesy.utily import _passarg # _MODS
38from math import asinh, atan, atan2, copysign, cosh, fabs, sin, sinh, sqrt
39try:
40 from math import exp2 as _exp2
41except ImportError: # Python 3.11-
43 def _exp2(x):
44 return pow(_2_0, x)
46__all__ = ()
47__version__ = '23.09.14'
49_TRIPS = 1024 # XXX 2 or 3?
52class AuxLat(AuxAngle):
53 '''Accurate conversion between I{Auxiliary} latitudes on an ellipsoid.
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}.
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 _mAL = 6 # 4, 6 or 8 aka Lmax
67 def __init__(self, a_earth=_WGS84, f=None, b=None, name=NN, **ALorder):
68 '''New L{AuxLat} instance on an ellipsoid or datum.
70 @arg a_earth: Equatorial radius, semi-axis (C{meter}) or an
71 ellipsoid or datum (L{Datum}, L{Ellipsoid},
72 L{Ellipsoid2} or L{a_f2Tuple}).
73 @kwarg f: Flattening: M{(a - b) / a} (C{float}, near zero for
74 spherical), ignored if B{C{a_earth}} is not scalar.
75 @kwarg b: Optional polar radius, semi-axis (C{meter}, same
76 units as B{C{a_earth}}), ignored if B{C{a_earth}}
77 is not scalar.
78 @kwarg ALorder: Optional keyword arguments B{C{ALorder}} to
79 set L{AuxLat}'s C{order}, see property
80 C{ALorder}.
81 @kwarg name: Optional name (C{str}).
82 '''
83 if a_earth is not _WGS84:
84 n = name or AuxLat.__name__
85 try:
86 if b is f is None:
87 E = _ellipsoidal_datum(a_earth, name=n).ellipsoid # XXX raiser=_earth_
88 elif isscalar(a_earth):
89 E = Ellipsoid(a_earth, f=f, b=b, name=_UNDER_(n))
90 else:
91 raise ValueError()
92 except Exception as x:
93 raise AuxError(a_earth=a_earth, f=f, b=b, cause=x)
94 self._E = E
96 if name:
97 self.name = name
98 if ALorder:
99 self.ALorder = _xkwds_get(ALorder, ALorder=AuxLat._mAL)
101 @Property_RO
102 def a(self):
103 '''Get the C{ellipsoid}'s equatorial radius (C{meter}, conventionally).
104 '''
105 return self.ellipsoid.a
107 equatoradius = a
109 @Property
110 def ALorder(self):
111 '''Get the I{AuxLat} order (C{int}, 4, 6 or 8).
112 '''
113 return self._mAL
115 @ALorder.setter # PYCHOK setter!
116 def ALorder(self, order):
117 '''Set the I{AuxLat} order (C{int}, 4, 6 or 8).
118 '''
119 m = _Xorder(_AR2Coeffs, AuxError, ALorder=order)
120 if self._mAL != m:
121 _update_all(self)
122 self._mAL = m
124 def _atanhee(self, tphi): # see Ellipsoid._es_atanh, .albers._atanhee
125 # atanh(e * sphi) = asinh(e' * sbeta)
126 f = self.f
127 s = _sn(self._fm1 * tphi) if f > 0 else _sn(tphi)
128 if f: # atanh(e * sphi) = asinh(e' * sbeta)
129 e = self._e
130 s = (atan(e * s) if f < 0 else asinh(self._e1 * s)) / e
131 return s
133 def Authalic(self, Phi, **diff_name):
134 '''Convert I{Geographic} to I{Aunthalic} latitude.
136 @arg Phi: Geographic latitude (L{AuxAngle}).
138 @return: Parametric latitude, C{Beta} (L{AuxAngle}).
139 '''
140 _xinstanceof(AuxAngle, Phi=Phi)
141 # assert Phi._AUX == Aux.PHI
142 tphi = fabs(Phi.tan)
143 if isfinite(tphi) and tphi and self.f:
144 y, x = Phi._yx_normalized
145 q = self._q
146 qv = self._qf(tphi)
147 Dq2 = self._Dq(tphi)
148 Dq2 *= (q + qv) / (fabs(y) + _1_0) # _Dq(-tphi)
149 Xi = AuxXi(copysign(qv, Phi.y), x * sqrt(Dq2),
150 name=_xkwds_get(diff_name, name=Phi.name))
152 if _xkwds_get(diff_name, diff=False):
153 if isnan(tphi):
154 d = self._e2m1_sq2
155 else:
156 c = self.Parametric(Phi)._x_normalized
157 d = _over(c, Xi._x_normalized)**3
158 d *= _over(c, x) * _over(_2_0, q)
159 Xi._diff = d
160 else:
161 Xi = AuxXi(*Phi._yx) # diff default
162 # assert Xi._AUX == Aux.XI
163 return Xi
165 def AuthalicRadius2(self, exact=False, f_max=_0_1):
166 '''Get the I{Authalic} radius I{squared}.
168 @kwarg exact: If C{True}, use the exact expression, otherwise
169 the I{Taylor} series.
170 @kwarg f_max: C{Flattening} not to exceed (C{float}).
172 @return: Authalic radius I{squared} (C{meter} I{squared}, same
173 units as the ellipsoid axes).
175 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)}
176 exceeds C{f_max}.
177 '''
178 f = self.f
179 if exact or not f:
180 c2 = self.ellipsoid.b2 * self._q # == ellipsoid.c2x * 2
181 elif fabs(f) < f_max:
182 # Using a * (a + b) / 2 as the multiplying factor leads to a rapidly
183 # converging series in n. Of course, using this series isn't really
184 # necessary, since the exact expression is simple to evaluate. However,
185 # we do it for consistency with RectifyingRadius and, presumably, the
186 # roundoff error is smaller compared to that for the exact expression.
187 m = self.ALorder
188 c2 = _MODS.karney._polynomial(self._n, _AR2Coeffs[m], 0, m)
189 c2 *= self.a * (self.a + self.b)
190 else:
191 raise AuxError(exact=exact, f=f, f_max=f_max)
192 return c2 * _0_5
194 @Property_RO
195 def b(self):
196 '''Get the C{ellipsoid}'s polar radius (C{meter}, conventionally).
197 '''
198 return self.ellipsoid.b # a * (_1_0 - f)
200 polaradius = b
202 def _coeffs(self, auxout, auxin):
203 # Get the polynomial coefficients as 4-, 6- or 8-tuple
204 aL = self.ALorder # aka Lmax
205 k = Aux._1d(auxout, auxin)
206 try:
207 cs = AuxLat._csc[aL][k]
208 except KeyError:
209 if auxout == auxin:
210 cs = _0_0s(aL) # not cached
211 else:
212 Cx = _CXcoeffs(aL)
213 try:
214 Cx = Cx[auxout][auxin]
215 except KeyError as x:
216 raise AuxError(auxout=auxout, auxin=auxin, cause=x)
218 d = x = n = self._n
219 if Aux.use_n2(auxin) and Aux.use_n2(auxout):
220 x = self._n2
222 def _m(m):
223 return m // 2
224 else:
225 def _m(m): # PYCHOK expected
226 return m
228 i = 0
229 cs = []
230 _p = _MODS.karney._polynomial
231 for r in _reverange(aL):
232 j = i + _m(r) + 1 # order m = j - i - 1
233 cs.append(_p(x, Cx, i, j) * d)
234 d *= n
235 i = j
236 # assert i == len(Cx) and len(cs) == aL
237 AuxLat._csc.setdefault(aL, {})[k] = tuple(cs)
238 return cs
240 def Conformal(self, Phi, **diff_name):
241 '''Convert I{Geographic} to I{Conformal} latitude.
243 @arg Phi: Geographic latitude (L{AuxAngle}).
245 @return: Conformal latitude, C{Chi} (L{AuxAngle}).
246 '''
247 _xinstanceof(AuxAngle, Phi=Phi)
248 # assert Phi._AUX == Aux.PHI
249 tphi = tchi = fabs(Phi.tan)
250 if isfinite(tphi) and tphi and self.f:
251 sig = sinh(self._e2 * self._atanhee(tphi))
252 scsig = _sc(sig)
253 scphi = _sc(tphi)
254 if self.f > 0:
255 # The general expression for tchi is
256 # tphi * scsig - sig * scphi
257 # This involves cancellation if f > 0, so change to
258 # (tphi - sig) * (tphi + sig) / (tphi * scsig + sig * scphi)
259 # To control overflow, write as (sigtphi = sig / tphi)
260 # (tphi - sig) * (1 + sigtphi) / (scsig + sigtphi * scphi)
261 sigtphi = sig / tphi
262 if sig < (tphi * _0_5):
263 t = tphi - sig
264 else:
265 def _asinh_2(x):
266 return asinh(x) * _0_5
267 # Still possibly dangerous cancellation in tphi - sig.
268 # Write tphi - sig = (1 - e) * Dg(1, e)
269 # Dg(x, y) = (g(x) - g(y)) / (x - y)
270 # g(x) = sinh(x * atanh(sphi * x))
271 # Note sinh(atanh(sphi)) = tphi
272 # Turn the crank on divided differences, substitute
273 # sphi = tphi / _sc(tphi)
274 # atanh(x) = asinh(x / sqrt(1 - x^2))
275 e = self._e
276 em1 = self._e2m1 / (_1_0 + e)
277 # assert em1 != 0
278 scb = self._scbeta(tphi)
279 scphib = scphi / scb # sec(phi) / sec(beta)
280 atphib = _asinh_2(tphi * e / scb) # atanh(e * sphi)
281 atphi = _asinh_2(tphi) # atanh(sphi)
282 t = _asinh_2(em1 * (tphi * scphib)) / em1
283 try:
284 Dg = Fsum(atphi, atphib, t, e * t)
285 except ValueError: # Fsum(NAN) exception
286 Dg = _sum((atphi, atphib, t, e * t))
287 e *= atphib
288 t = atphi - e
289 if t: # sinh(0) == 0
290 Dg *= sinh(t) / t * cosh(atphi + e) * em1
291 t = float(Dg) # tphi - sig
292 tchi = _over(t * (_1_0 + sigtphi),
293 scsig + scphi * sigtphi) # if t else _0_0
294 else:
295 tchi = tphi * scsig - sig * scphi
297 n = _xkwds_get(diff_name, name=Phi.name)
298 Chi = AuxChi(tchi, name=n).copyquadrant(Phi)
300 if _xkwds_get(diff_name, diff=False):
301 if isinf(tphi): # PYCHOK np cover
302 d = self._conformal_diff
303 else:
304 d = self.Parametric(Phi)._x_normalized
305 if d:
306 d = _over(d, Chi._x_normalized) * \
307 _over(d, Phi._x_normalized) * self._e2m1
308 Chi._diff = d
309 # assrt Chi._AUX == Aux.CHI
310 return Chi
312 @Property_RO
313 def _conformal_diff(self): # PYCHOK no cover
314 '''(INTERNAL) Constant I{Conformal} diff.
315 '''
316 e = self._e
317 if self.f > 0:
318 ss = sinh(asinh(self._e1) * e)
319 d = _over(_1_0, _sc(ss) + ss)
320 elif e:
321 ss = sinh(-atan(e) * e)
322 d = _sc(ss) - ss
323 else:
324 d = _1_0
325 return d
327 def convert(self, auxout, Zeta_d, exact=False):
328 # Convert I{Auxiliary} or I{scalar} latitude
329 Z = d = Zeta_d
330 if isinstance(Z, AuxAngle):
331 A, auxin = _AuxClass(auxout), Z._AUX
332 if auxin == auxout:
333 pass
334 elif not (isfinite(Z.tan) and Z.tan): # XXX
335 Z = A(*Z._yx, aux=auxout, name=Z.name)
336 elif exact:
337 p = Aux.power(auxout, auxin)
338 if p is None:
339 P = self._fromAux(Z) # Phi
340 Z = self._toAux(auxout, P)
341 Z._iter = P.iteration
342 else:
343 y, x = Z._yx
344 if p:
345 y *= pow(self._fm1, p)
346 Z = A(y, x, aux=auxout, name=Z.name)
347 else:
348 cs = self._coeffs(auxout, auxin)
349 yx = Z._yx_normalized
350 Z = A(*yx, aux=auxout, name=Z.name)
351 # assert Z._yx == yx
352 r = _Clenshaw(True, Z, cs, self.ALorder)
353 Z += AuxAngle.fromRadians(r, aux=auxout)
354 # assert Z._AUX == auxout
355 return Z
357 elif isinstance(d, Degrees) or isscalar(d):
358 Z = AuxPhi.fromDegrees(d)
359 d = round((d - Z.toDegrees) / _360_0) * _360_0
360 d += self.convert(auxout, Z, exact).toDegrees
361 return Degrees(d, name=Aux.Greek(auxout))
363 raise AuxError(auxout=auxout, Zeta_d=Zeta_d, exact=exact)
365 def _Dq(self, tphi):
366 # I{Divided Difference} of (q(1) - q(sphi)) / (1 - sphi).
367 sphi = _sn(tphi)
368 if tphi > 0:
369 scphi = _sc(tphi)
370 d = _1_0 / (scphi**2 * (_1_0 + sphi)) # XXX - sphi
371 if d:
372 # General expression for _Dq(1, sphi) is
373 # atanh(e * d / (1 - e2 * sphi)) / (e * d) +
374 # (1 + e2 * sphi) / ((1 - e2 * sphi * sphi) * e2m1)
375 # with atanh(e * d / (1 - e2 * sphi)) =
376 # atanh(e * d * scphi / (scphi - e2 * tphi))
377 e2m1, ed = self._e2m1, (self._e * d)
378 if e2m1 and ed:
379 e2 = self._e2
380 if e2 > 0: # assert self.f > 0
381 scb = self._scbeta(tphi)
382 q = scphib = scphi / scb
383 q *= (scphi + tphi * e2) / (e2m1 * scb)
384 q += asinh(self._e1 * d * scphib) / ed
385 else:
386 s2 = sphi * e2
387 q = (_1_0 + s2) / ((_1_0 - sphi * s2) * e2m1)
388 q += (atan2(ed, _1_0 - s2) / ed) if e2 < 0 else _1_0
389 else: # PYCHOK no cover
390 q = INF
391 else: # PYCHOK no cover
392 q = self._2_e2m12
393 else: # not reached, open-coded in .Authalic
394 q = _over(self._q - self._qf(tphi), _1_0 - sphi)
395 return q
397 @Property_RO
398 def _e(self): # unsigned, (1st) eccentricity
399 return self.ellipsoid.e # sqrt(fabs(self._e2))
401 @Property_RO
402 def _e1(self):
403 return sqrt(fabs(self._e12))
405 @Property_RO
406 def _e12(self):
407 return _over(self._e2, _1_0 - self._e2)
409 @Property_RO
410 def _e12p1(self):
411 return _1_0 / self._e2m1
413 @Property_RO
414 def _e2(self): # signed, (1st) eccentricity squared
415 return self.ellipsoid.e2
417 @Property_RO
418 def _e2m1(self): # 1 less 1st eccentricity squared
419 return self.ellipsoid.e21 # == ._fm1**2
421 @Property_RO
422 def _e2m1_sq2(self):
423 return self._e2m1 * sqrt(self._q * _0_5)
425 @Property_RO
426 def _2_e2m12(self):
427 return _2_0 / self._e2m1**2
429 @Property_RO
430 def _Ef_fRG_a2b2_PI_4(self):
431 E = self.ellipsoid
432 return _Ef.fRG(E.a2, E.b2) / PI_4
434 @Property_RO
435 def ellipsoid(self):
436 '''Get the ellipsoid (L{Ellipsoid}).
437 '''
438 return self._E
440 @Property_RO
441 def f(self):
442 '''Get the C{ellipsoid}'s flattening (C{scalar}).
443 '''
444 return self.ellipsoid.f
446 flattening = f
448 @Property_RO
449 def _fm1(self): # 1 - flattening
450 return self.ellipsoid.f1
452 def _fromAux(self, Zeta, **name):
453 '''Convert I{Auxiliary} to I{Geographic} latitude.
455 @arg Zeta: Auxiliary latitude (L{AuxAngle}).
457 @return: Geographic latitude, I{Phi} (L{AuxAngle}).
458 '''
459 _xinstanceof(AuxAngle, Zeta=Zeta)
460 aux = Zeta._AUX
461 n = _xkwds_get(name, name=Zeta.name)
462 f = self._fromAuxCase.get(aux, None)
463 if f is None:
464 Phi = AuxPhi(NAN, name=n)
465 elif callable(f):
466 t = self._fm1
467 t *= f(t)
468 Phi = _Newton(t, Zeta, self._toZeta(aux), name=n)
469 else: # assert isscalar(f)
470 y, x = Zeta._yx
471 Phi = AuxPhi(y / f, x, name=n)
472 # assert Phi._AUX == Aux.PHI
473 return Phi
475 @Property_RO
476 def _fromAuxCase(self):
477 '''(INTERNAL) switch(auxin): ...
478 '''
479 return {Aux.AUTHALIC: cbrt,
480 Aux.CONFORMAL: _MODS.utily._passarg,
481 Aux.GEOCENTRIC: self._e2m1,
482 Aux.GEOGRAPHIC: _1_0,
483 Aux.PARAMETRIC: self._fm1,
484 Aux.RECTIFYING: sqrt}
486 def Geocentric(self, Phi, **diff_name):
487 '''Convert I{Geographic} to I{Geocentric} latitude.
489 @arg Phi: Geographic latitude (L{AuxAngle}).
491 @return: Geocentric latitude, C{Phi} (L{AuxAngle}).
492 '''
493 _xinstanceof(AuxAngle, Phi=Phi)
494 # assert Phi._AUX == Aux.PHI
495 Theta = AuxTheta(Phi.y * self._e2m1, Phi.x,
496 name=_xkwds_get(diff_name, name=Phi.name))
497 if _xkwds_get(diff_name, diff=False):
498 Theta._diff = self._e2m1
499 return Theta
501 def Geodetic(self, Phi, **diff_name): # PYCHOK no cover
502 '''Convert I{Geographic} to I{Geodetic} latitude.
504 @arg Phi: Geographic latitude (L{AuxAngle}).
506 @return: Geodetic latitude, C{Phi} (L{AuxAngle}).
507 '''
508 _xinstanceof(AuxAngle, Phi=Phi)
509 # assert Phi._AUX == Aux.PHI
510 return AuxPhi(Phi, name=_xkwds_get(diff_name, name=Phi.name))
512 @Property_RO
513 def _n(self): # 3rd flattening
514 return self.ellipsoid.n
516 @Property_RO
517 def _n2(self):
518 return self._n**2
520 def Parametric(self, Phi, **diff_name):
521 '''Convert I{Geographic} to I{Parametric} latitude.
523 @arg Phi: Geographic latitude (L{AuxAngle}).
525 @return: Parametric latitude, C{Beta} (L{AuxAngle}).
526 '''
527 _xinstanceof(AuxAngle, Phi=Phi)
528 # assert Phi._AUX == Aux.PHI
529 Beta = AuxBeta(Phi.y * self._fm1, Phi.x,
530 name=_xkwds_get(diff_name, name=Phi.name))
531 if _xkwds_get(diff_name, diff=False):
532 Beta._diff = self._fm1
533 return Beta
535 Reduced = Parametric
537 @Property_RO
538 def _q(self): # constant _q
539 q, f = self._e12p1, self.f
540 if f:
541 e = self._e
542 q += (asinh(self._e1) if f > 0 else atan(e)) / e
543 else:
544 q += _1_0
545 return q
547 def _qf(self, tphi):
548 # function _q: atanh(e * sphi) / e + sphi / (1 - (e * sphi)^2)
549 scb = self._scbeta(tphi)
550 return self._atanhee(tphi) + (tphi / scb) * (_sc(tphi) / scb)
552 def _qIntegrand(self, beta):
553 # pbeta(beta) = integrate(q(beta), beta), with beta in radians
554 # q(beta) = (1-f) * (sin(xi) - sin(chi)) / cos(phi)
555 # = (1-f) * (cos(chi) - cos(xi)) / cos(phi) *
556 # (cos(xi) + cos(chi)) / (sin(xi) + sin(chi))
557 # Fit q(beta)/cos(beta) with Fourier transform
558 # q(beta)/cos(beta) = sum(c[k] * sin((2*k+1)*beta), k, 0, K-1)
559 # then the integral is
560 # pbeta = sum(d[k] * cos((2*k+2)*beta), k, 0, K-1)
561 # where
562 # d[k] = -1/(4*(k+1)) * (c[k] + c[k+1]) for k in 0..K-2
563 # d[K-1] = -1/(4*K) * c[K-1]
564 Beta = AuxBeta.fromRadians(beta)
565 if Beta.x: # and self._fm1:
566 Ax, _cv = Aux, self.convert
567 Phi = _cv(Ax.PHI, Beta, exact=True)
568 schi, cchi = _cv(Ax.CHI, Phi, exact=True)._yx_normalized
569 sxi, cxi = _cv(Ax.XI, Phi, exact=True)._yx_normalized
570 r = (sxi - schi) if fabs(schi) < fabs(cchi) else \
571 _over(_2cos2x(cchi, cxi), (sxi + schi) * _2_0)
572 r *= _over(self._fm1, Phi._x_normalized * Beta._x_normalized)
573 else: # beta == PI_2, PI3_2, ...
574 r = _0_0 # XXX 0 avoids NAN summation exceptions
575 return r
577 def Rectifying(self, Phi, **diff_name):
578 '''Convert I{Geographic} to I{Rectifying} latitude.
580 @arg Phi: Geographic latitude (L{AuxAngle}).
582 @return: Rectifying latitude, C{Mu} (L{AuxAngle}).
583 '''
584 Beta = self.Parametric(Phi)
585 # assert Beta._AUX == Aux.BETA
586 sb, cb = map(fabs, Beta._yx_normalized)
587 a, ka, ka1 = _1_0, self._e2, self._e2m1
588 b, kb, kb1 = self._fm1, -self._e12, self._e12p1
589 if self.f < 0:
590 a, b = b, a
591 ka, kb = kb, ka
592 ka1, kb1 = kb1, ka1
593 sb, cb = cb, sb
594 # now a, b = larger/smaller semiaxis
595 # Beta measured from larger semiaxis
596 # kb, ka = modulus-squared for distance from Beta = 0, pi/2
597 # NB kb <= 0; 0 <= ka <= 1
598 # sa = b*E(Beta, sqrt(kb))
599 # sb = a*E(Beta',sqrt(ka))
600 # 1 - ka * (1 - sb2) = 1 - ka + ka*sb2
601 sb2 = sb**2
602 cb2 = cb**2
603 da2 = ka1 + ka * sb2
604 db2 = _1_0 - kb * sb2
605 # DLMF Eq. 19.25.9
606 my = b * sb * _Ef._RFRD(cb2, db2, _1_0, kb * sb2)
607 # DLMF Eq. 19.25.10 with complementary angles
608 mx = a * cb * (_Ef.fRF(sb2, da2, _1_0) * ka1 +
609 ka * cb2 * _Ef.fRD(sb2, _1_0, da2, _3_0) * ka1 +
610 ka * sb / sqrt(da2))
611 # my + mx = 2*_Ef.fRG(a*a, b*b) = a*E(e) = b*E(i*e')
612 # mr = a*E(e)*(2/pi) = b*E(i*e')*(2/pi)
613 if self.f < 0:
614 a, b = b, a
615 my, mx = mx, my
616 mr = (my + mx) / PI_2
617 if mr:
618 my = sin(my / mr)
619 mx = sin(mx / mr) # XXX zero?
620 else: # zero Mu
621 my, mx = _0_0, _1_0
622 n = _xkwds_get(diff_name, name=Phi.name)
623 Mu = AuxMu(my, mx, # normalized
624 name=n).copyquadrant(Phi)
626 if _xkwds_get(diff_name, diff=False):
627 d, x = _0_0, Beta._x_normalized
628 if x and mr:
629 if Mu.x and Phi.x and not isinf(Phi.tan):
630 d = b / mr * (x / Mu.x)**2 \
631 * (x / Phi._x_normalized)
632 else:
633 d = mr / a
634 Mu._diff = self._fm1 * d
635 return Mu
637 def RectifyingRadius(self, exact=False):
638 '''Get the I{Rectifying} radius.
640 @arg exact: If C{True}, use the exact expression,
641 otherwise the I{Taylor} series.
643 @return: Rectifying radius (L{Meter}, same units
644 as the ellipsoid axes).
645 '''
646 r = self._Ef_fRG_a2b2_PI_4 if exact else self._RectifyingR
647 return Meter(r, name=self.RectifyingRadius.__name__)
649 @Property_RO
650 def _RectifyingR(self):
651 m = self.ALorder
652 d = _MODS.karney._polynomial(self._n2, _RRCoeffs[m], 0, m // 2)
653 return d * (self.a + self.b) * _0_5
655 def _scbeta(self, tphi):
656 return _sc(self._fm1 * tphi)
658 def _toAux(self, auxout, Phi, **diff_name):
659 '''Convert I{Geographic} to I{Auxiliary} latitude.
661 @arg auxout: I{Auxiliary} kind (C{Aux.KIND}).
662 @arg Phi: Geographic latitude (L{AuxLat}).
664 @return: Auxiliary latitude, I{Eta} (L{AuxLat}).
665 '''
666 _xinstanceof(AuxAngle, Phi=Phi)
667 # assert Phi._AUX == Aux.PHI
668 n = _xkwds_get(diff_name, name=Phi.name)
669 m = _toAuxCase.get(auxout, None)
670 if m: # callable
671 A = m(self, Phi, **_xkwds(diff_name, name=n))
672 elif auxout == Aux.GEODETIC: # == GEOGRAPHIC
673 A = AuxPhi(Phi, name=n)
674 else: # auxout?
675 A = AuxPhi(NAN, name=n)
676 # assert A._AUX == auxout
677 return A
679 def _toZeta(self, zetaux):
680 '''Return a (lean) function to create C{AuxPhi(tphi)} and
681 convert that into C{AuxAngle} of (fixed) kind C{zetaux}
682 for use only inside the C{_Newton} loop.
683 '''
684 class _AuxPhy(AuxPhi):
685 # lean C{AuxPhi} instance.
686 # _diff = _1_0
687 # _x = _1_0
689 def __init__(self, tphi): # PYCHOK signature
690 self._y = tphi
692 m = _toAuxCase.get(zetaux, None)
693 if m: # callable
695 def _toZeta(tphi):
696 return m(self, _AuxPhy(tphi), diff=True)
698 elif zetaux == Aux.GEODETIC: # GEOGRAPHIC
699 _toZeta = _AuxPhy
701 else: # zetaux?
703 def _toZeta(unused): # PYCHOK expected
704 return _AuxPhy(NAN)
706 return _toZeta
709# switch(auxout): ...
710_toAuxCase = {Aux.AUTHALIC: AuxLat.Authalic,
711 Aux.CONFORMAL: AuxLat.Conformal,
712 Aux.GEOCENTRIC: AuxLat.Geocentric,
713 Aux.PARAMETRIC: AuxLat.Parametric,
714 Aux.RECTIFYING: AuxLat.Rectifying}
717def _Clenshaw(sinp, Zeta, cs, K):
718 sz, cz = Zeta._yx # isnormal
719 # Evaluate sum(c[k] * sin((2*k+2) * Zeta)) if sinp else
720 # sum(c[k] * cos((2*k+2) * Zeta))
721 x = _2cos2x(cz, sz) # 2 * cos(2*Zeta)
722 if isfinite(x):
723 U0, U1 = Fsum(), Fsum()
724 # assert len(cs) == K
725 for r in _reverange(K):
726 U1 -= U0 * x + cs[r]
727 U1, U0 = U0, -U1
728 # u0*f0(Zeta) - u1*fm1(Zeta)
729 # f0 = sin(2*Zeta) if sinp else cos(2*Zeta)
730 # fm1 = 0 if sinp else 1
731 if sinp:
732 U0 *= sz * cz * _2_0
733 else:
734 U0 *= x * _0_5
735 U0 -= U1
736 x = float(U0)
737 return x
740def _CXcoeffs(aL): # PYCHOK in .auxilats.__main__
741 '''(INTERNAL) Get the C{CX_4}, C{_6} or C{_8} coefficients.
742 '''
743 try: # from pygeodesy.auxilats._CX_x import _coeffs_x as _coeffs
744 _CX_x = _DOT_(_MODS.auxilats.__name__, _UNDER_('_CX', aL))
745 _coeffs = _MODS.getattr(_CX_x, _UNDER_('_coeffs', aL))
746 except (AttributeError, ImportError, KeyError, TypeError) as x:
747 raise AuxError(ALorder=aL, cause=x)
748 # assert _coeffs.ALorder == aL
749 # assert _coeffs.n == Aux.len(aL)
750 return _coeffs
753def _Newton(tphi, Zeta, _toZeta, **name):
754 # Newton's method fro AuxLat._fromAux
755 try:
756 _lg2 = _log2
757 _abs = fabs
758 tz = _abs(Zeta.tan)
759 tphi = tz / tphi # **)
760 ltz = _lg2(tz) # **)
761 ltphi = _lg2(tphi) # **)
762 ltmin = min(ltphi, MIN_EXP)
763 ltmax = max(ltphi, MAX_EXP)
764# auxin = Zeta._AUX
765 s, n, __2 = 0, 3, _0_5 # n = i + 2
766 _TOL, _xp2 = _EPSqrt, _exp2
767 for i in range(1, _TRIPS): # up to 1 Ki!
768 # _toAux(auxin, AuxPhi(tphi), diff=True)
769 Z = _toZeta(tphi)
770 # assert Z._AUX == auxin
771 t, s_ = Z.tan, s
772 if t > tz:
773 ltmax, s = ltphi, +1
774 elif t < tz:
775 ltmin, s = ltphi, -1
776 else:
777 break
778 # derivative dtan(Z)/dtan(Phi)
779 # to dlog(tan(Z))/dlog(tan(Phi))
780 d = (ltz - _lg2(t)) * t # **)
781 if d:
782 d = d / (Z.diff * tphi) # **)
783 ltphi += d
784 tphi = _xp2(ltphi)
785 if _abs(d) < _TOL:
786 i += 1
787 # _toAux(auxin, AuxPhi(tphi), diff=True)
788 Z = _toZeta(tphi)
789 tphi -= _over(Z.tan - tz, Z.diff)
790 break
791 if (i > n and (s * s_) < 0) or not ltmin < ltphi < ltmax:
792 s, n = 0, (i + 2)
793 ltphi = (ltmin + ltmax) * __2
794 tphi = _xp2(ltphi)
795 else:
796 i = _TRIPS
797 Phi = AuxPhi(tphi, **name).copyquadrant(Zeta)
798 Phi._iter = i
799 except (ValueError, ZeroDivisionError): # **) zero t, tphi, tz or Z.diff
800 Phi = AuxPhi(Zeta, **name) # diff as-as
801 Phi._iter = 0
802 # assert Phi._AUX == Aux.PHI
803 return Phi
806_f, _u = float, _Ufloats()
807_1__f3 = -1 / _f(3) # XXX +1 / _f(3)
808_AR2Coeffs = {4: _u(4 / _f(315), 4 / _f(105), 4 / _f(15), _1__f3),
809 6: _u(4 / _f(1287), 4 / _f(693), 4 / _f(315), 4 / _f(105),
810 4 / _f(15), _1__f3),
811 8: _u(4 / _f(3315), 4 / _f(2145), 4 / _f(1287), 4 / _f(693),
812 4 / _f(315), 4 / _f(105), 4 / _f(15), _1__f3)}
813_RRCoeffs = {4: _u(1 / _f(64), _0_25),
814 6: _u(1 / _f(256), 1 / _f(64), _0_25),
815 8: _u(25 / _f(16384), 1 / _f(256), 1 / _f(64), _0_25)} # PYCHOK used!
816del _f, _u, _Ufloats, _1__f3
817# assert set(_AR2Coeffs.keys()) == set(_RRCoeffs.keys())
819__all__ += _ALL_DOCS(AuxLat)
821# **) MIT License
822#
823# Copyright (C) 2023-2023 -- mrJean1 at Gmail -- All Rights Reserved.
824#
825# Permission is hereby granted, free of charge, to any person obtaining a
826# copy of this software and associated documentation files (the "Software"),
827# to deal in the Software without restriction, including without limitation
828# the rights to use, copy, modify, merge, publish, distribute, sublicense,
829# and/or sell copies of the Software, and to permit persons to whom the
830# Software is furnished to do so, subject to the following conditions:
831#
832# The above copyright notice and this permission notice shall be included
833# in all copies or substantial portions of the Software.
834#
835# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
836# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
837# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
838# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
839# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
840# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
841# OTHER DEALINGS IN THE SOFTWARE.