Coverage for pygeodesy/geodesicx/gx.py: 93%
524 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-25 12:04 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-25 12:04 -0400
2# -*- coding: utf-8 -*-
4u'''A pure Python version of I{Karney}'s C++ class U{GeodesicExact
5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}.
7Class L{GeodesicExact} follows the naming, methods and return values
8of class C{Geodesic} from I{Karney}'s Python U{geographiclib
9<https://GitHub.com/geographiclib/geographiclib-python>}.
11Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023)
12and licensed under the MIT/X11 License. For more information, see the
13U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
14'''
15# make sure int/int division yields float quotient
16from __future__ import division as _; del _ # PYCHOK semicolon
18# A copy of comments from Karney's C{GeodesicExact.cpp}:
19#
20# This is a reformulation of the geodesic problem. The
21# notation is as follows:
22# - at a general point (no suffix or 1 or 2 as suffix)
23# - phi = latitude
24# - lambda = longitude
25# - beta = latitude on auxiliary sphere
26# - omega = longitude on auxiliary sphere
27# - alpha = azimuth of great circle
28# - sigma = arc length along great circle
29# - s = distance
30# - tau = scaled distance (= sigma at multiples of PI/2)
31# - at northwards equator crossing
32# - beta = phi = 0
33# - omega = lambda = 0
34# - alpha = alpha0
35# - sigma = s = 0
36# - a 12 suffix means a difference, e.g., s12 = s2 - s1.
37# - s and c prefixes mean sin and cos
39from pygeodesy.basics import _copysign, _xinstanceof, _xor, unsigned0
40from pygeodesy.constants import EPS, EPS0, EPS02, MANT_DIG, NAN, PI, _EPSqrt, \
41 _SQRT2_2, isnan, _0_0, _0_001, _0_01, _0_1, _0_5, \
42 _1_0, _N_1_0, _1_75, _2_0, _N_2_0, _2__PI, _3_0, \
43 _4_0, _6_0, _8_0, _16_0, _90_0, _180_0, _1000_0
44from pygeodesy.datums import _earth_datum, _WGS84, _EWGS84
45# from pygeodesy.ellipsoids import _EWGS84 # from .datums
46from pygeodesy.errors import GeodesicError, _xkwds_pop2
47from pygeodesy.fmath import hypot as _hypot
48from pygeodesy.fsums import fsumf_, fsum1f_
49from pygeodesy.geodesicx.gxbases import _cosSeries, _GeodesicBase, \
50 _sincos12, _sin1cos2, _xnC4
51from pygeodesy.geodesicx.gxline import _GeodesicLineExact, _TINY, _update_glXs
52from pygeodesy.interns import NN, _COMMASPACE_, _DOT_, _UNDER_
53from pygeodesy.karney import GDict, _around, _atan2d, Caps, _cbrt, _diff182, \
54 _fix90, _K_2_0, _norm2, _norm180, _polynomial, \
55 _signBit, _sincos2, _sincos2d, _sincos2de, _unsigned2
56from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS
57from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple
58from pygeodesy.props import deprecated_Property, Property, Property_RO, property_RO
59from pygeodesy.streprs import Fmt, pairs
60from pygeodesy.utily import atan2d as _atan2d_reverse, _unrollon, _Wrap, wrap360
62from math import atan2, copysign, cos, degrees, fabs, radians, sqrt
64__all__ = ()
65__version__ = '24.05.20'
67_MAXIT1 = 20
68_MAXIT2 = 10 + _MAXIT1 + MANT_DIG # MANT_DIG == C++ digits
70# increased multiplier in defn of _TOL1 from 100 to 200 to fix Inverse
71# case 52.784459512564 0 -52.784459512563990912 179.634407464943777557
72# which otherwise failed for Visual Studio 10 (Release and Debug)
73_TOL0 = EPS
74_TOL1 = _TOL0 * -200 # negative
75_TOL2 = _EPSqrt # == sqrt(_TOL0)
76_TOL3 = _TOL2 * _0_1
77_TOLb = _TOL2 * _TOL0 # Check on bisection interval
78_THR1 = _TOL2 * _1000_0 + _1_0
80_TINY3 = _TINY * _3_0
81_TOL08 = _TOL0 * _8_0
82_TOL016 = _TOL0 * _16_0
85def _atan12(*sincos12, **sineg0):
86 '''(INTERNAL) Return C{ang12} in C{radians}.
87 '''
88 return atan2(*_sincos12(*sincos12, **sineg0))
91def _eTOL2(f):
92 # Using the auxiliary sphere solution with dnm computed at
93 # (bet1 + bet2) / 2, the relative error in the azimuth
94 # consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
95 # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000.
97 # For a given f and sig12, the max error occurs for lines
98 # near the pole. If the old rule for computing dnm = (dn1
99 # + dn2)/2 is used, then the error increases by a factor of
100 # 2.) Setting this equal to epsilon gives sig12 = etol2.
102 # Here 0.1 is a safety factor (error decreased by 100) and
103 # max(0.001, abs(f)) stops etol2 getting too large in the
104 # nearly spherical case.
105 t = min(_1_0, _1_0 - f * _0_5) * max(_0_001, fabs(f)) * _0_5
106 return _TOL3 / (sqrt(t) if t > EPS02 else EPS0)
109class _PDict(GDict):
110 '''(INTERNAL) Parameters passed around in C{._GDictInverse} and
111 optionally returned when C{GeodesicExact.debug} is C{True}.
112 '''
113 def set_sigs(self, ssig1, csig1, ssig2, csig2):
114 '''Update the C{sig1} and C{sig2} parameters.
115 '''
116 self.set_(ssig1=ssig1, csig1=csig1, sncndn1=(ssig1, csig1, self.dn1), # PYCHOK dn1
117 ssig2=ssig2, csig2=csig2, sncndn2=(ssig2, csig2, self.dn2)) # PYCHOK dn2
119 def toGDict(self): # PYCHOK no cover
120 '''Return as C{GDict} without attrs C{sncndn1} and C{sncndn2}.
121 '''
122 def _rest(sncndn1=None, sncndn2=None, **rest): # PYCHOK sncndn* not used
123 return GDict(rest)
125 return _rest(**self)
128class GeodesicExact(_GeodesicBase):
129 '''A pure Python version of I{Karney}'s C++ class U{GeodesicExact
130 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>},
131 modeled after I{Karney}'s Python class U{geodesic.Geodesic<https://GitHub.com/
132 geographiclib/geographiclib-python>}.
133 '''
134 _datum = _WGS84
135 _nC4 = 30 # default C4order
137 def __init__(self, a_ellipsoid=_EWGS84, f=None, C4order=None, **name_C4Order): # for backward compatibility
138 '''New L{GeodesicExact} instance.
140 @arg a_ellipsoid: An ellipsoid (L{Ellipsoid}) or datum (L{Datum}) or
141 the equatorial radius of the ellipsoid (C{scalar},
142 conventionally in C{meter}), see B{C{f}}.
143 @arg f: The flattening of the ellipsoid (C{scalar}) if B{C{a_ellipsoid}}
144 is specified as C{scalar}.
145 @kwarg C4order: Optional series expansion order (C{int}), see property
146 L{C4order}, default C{30}.
147 @kwarg name_C4Order: Optional C{B{name}=NN} (C{str}) and the DEPRECATED
148 keyword argument C{C4Order}, use B{C{C4order}} instead.
150 @raise GeodesicError: Invalid B{C{C4order}}.
151 '''
152 if name_C4Order:
153 C4Order, name = _xkwds_pop2(name_C4Order, C4Order=C4order)
154 if C4Order: # for backward compatibility
155 self.C4order = C4Order
156 if name:
157 self.name = name
158 else:
159 name = {} # name_C4Order
161 _earth_datum(self, a_ellipsoid, f=f, **name)
162 if C4order: # XXX private copy, always?
163 self.C4order = C4order
165 @Property_RO
166 def a(self):
167 '''Get the I{equatorial} radius, semi-axis (C{meter}).
168 '''
169 return self.ellipsoid.a
171 def ArcDirect(self, lat1, lon1, azi1, a12, outmask=Caps.STANDARD):
172 '''Solve the I{Direct} geodesic problem in terms of (spherical) arc length.
174 @arg lat1: Latitude of the first point (C{degrees}).
175 @arg lon1: Longitude of the first point (C{degrees}).
176 @arg azi1: Azimuth at the first point (compass C{degrees}).
177 @arg a12: Arc length between the points (C{degrees}), can be negative.
178 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying
179 the quantities to be returned.
181 @return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2,
182 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1},
183 C{lon1}, C{azi1} and arc length C{a12} always included.
185 @see: C++ U{GeodesicExact.ArcDirect
186 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}
187 and Python U{Geodesic.ArcDirect<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
188 '''
189 return self._GDictDirect(lat1, lon1, azi1, True, a12, outmask)
191 def ArcDirectLine(self, lat1, lon1, azi1, a12, caps=Caps.ALL, **name):
192 '''Define a L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as arc length.
194 @arg lat1: Latitude of the first point (C{degrees}).
195 @arg lon1: Longitude of the first point (C{degrees}).
196 @arg azi1: Azimuth at the first point (compass C{degrees}).
197 @arg a12: Arc length between the points (C{degrees}), can be negative.
198 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying
199 the capabilities the L{GeodesicLineExact} instance
200 should possess, i.e., which quantities can be
201 returned by calls to L{GeodesicLineExact.Position}
202 and L{GeodesicLineExact.ArcPosition}.
203 @kwarg name: Optional C{B{name}=NN} (C{str}).
205 @return: A L{GeodesicLineExact} instance.
207 @note: The third point of the L{GeodesicLineExact} is set to correspond
208 to the second point of the I{Inverse} geodesic problem.
210 @note: Latitude B{C{lat1}} should in the range C{[-90, +90]}.
212 @see: C++ U{GeodesicExact.ArcDirectLine
213 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and
214 Python U{Geodesic.ArcDirectLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
215 '''
216 return self._GenDirectLine(lat1, lon1, azi1, True, a12, caps, **name)
218 def Area(self, polyline=False, **name):
219 '''Set up a L{GeodesicAreaExact} to compute area and
220 perimeter of a polygon.
222 @kwarg polyline: If C{True} perimeter only, otherwise
223 area and perimeter (C{bool}).
224 @kwarg name: Optional C{B{name}=NN} (C{str}).
226 @return: A L{GeodesicAreaExact} instance.
228 @note: The B{C{debug}} setting is passed as C{verbose}
229 to the returned L{GeodesicAreaExact} instance.
230 '''
231 gaX = _MODS.geodesicx.GeodesicAreaExact(self, polyline=polyline,
232 name=self._name__(name))
233 if self.debug:
234 gaX.verbose = True
235 return gaX
237 @Property_RO
238 def b(self):
239 '''Get the ellipsoid's I{polar} radius, semi-axis (C{meter}).
240 '''
241 return self.ellipsoid.b
243 @Property_RO
244 def c2x(self):
245 '''Get the ellipsoid's I{authalic} earth radius I{squared} (C{meter} I{squared}).
246 '''
247 # The Geodesic class substitutes atanh(sqrt(e2)) for asinh(sqrt(ep2))
248 # in the definition of _c2. The latter is more accurate for very
249 # oblate ellipsoids (which the Geodesic class does not handle). Of
250 # course, the area calculation in GeodesicExact is still based on a
251 # series and only holds for moderately oblate (or prolate) ellipsoids.
252 return self.ellipsoid.c2x
254 c2 = c2x # in this particular case
256 def C4f(self, eps):
257 '''Evaluate the C{C4x} coefficients for B{C{eps}}.
259 @arg eps: Polynomial factor (C{float}).
261 @return: C{C4order}-Tuple of C{C4x(B{eps})} coefficients.
262 '''
263 def _c4(nC4, C4x):
264 i, x, e = 0, _1_0, eps
265 _p = _polynomial
266 for r in range(nC4, 0, -1):
267 j = i + r
268 yield _p(e, C4x, i, j) * x
269 x *= e
270 i = j
271 # assert i == (nC4 * (nC4 + 1)) // 2
273 return tuple(_c4(self._nC4, self._C4x))
275 def _C4f_k2(self, k2): # in ._GDictInverse and gxline._GeodesicLineExact._C4a
276 '''(INTERNAL) Compute C{eps} from B{C{k2}} and invoke C{C4f}.
277 '''
278 return self.C4f(k2 / fsumf_(_2_0, sqrt(k2 + _1_0) * _2_0, k2))
280 @Property
281 def C4order(self):
282 '''Get the series expansion order (C{int}, 24, 27 or 30).
283 '''
284 return self._nC4
286 @C4order.setter # PYCHOK .setter!
287 def C4order(self, order):
288 '''Set the series expansion order (C{int}, 24, 27 or 30).
290 @raise GeodesicError: Invalid B{C{order}}.
291 '''
292 _xnC4(C4order=order)
293 if self._nC4 != order:
294 GeodesicExact._C4x._update(self)
295 _update_glXs(self) # zap cached _GeodesicLineExact attrs _B41, _C4a
296 self._nC4 = order
298 @deprecated_Property
299 def C4Order(self):
300 '''DEPRECATED, use property C{C4order}.
301 '''
302 return self.C4order
304 @C4Order.setter # PYCHOK .setter!
305 def C4Order(self, order):
306 '''DEPRECATED, use property C{C4order}.
307 '''
308 _xnC4(C4Order=order)
309 self.C4order = order
311 @Property_RO
312 def _C4x(self):
313 '''Get this ellipsoid's C{C4} coefficients, I{cached} tuple.
315 @see: Property L{C4order}.
316 '''
317 # see C4coeff() in GeographicLib.src.GeodesicExactC4.cpp
318 def _C4(nC4):
319 i, n, cs = 0, self.n, _C4coeffs(nC4)
320 _p = _polynomial
321 for r in range(nC4 + 1, 1, -1):
322 for j in range(1, r):
323 j = j + i # (j - i - 1) order of polynomial
324 yield _p(n, cs, i, j) / cs[j]
325 i = j + 1
326 # assert i == (nC4 * (nC4 + 1) * (nC4 + 5)) // 6
328 return tuple(_C4(self._nC4)) # 3rd flattening
330 @property_RO
331 def datum(self):
332 '''Get the datum (C{Datum}).
333 '''
334 return self._datum
336 def Direct(self, lat1, lon1, azi1, s12=0, outmask=Caps.STANDARD):
337 '''Solve the I{Direct} geodesic problem
339 @arg lat1: Latitude of the first point (C{degrees}).
340 @arg lon1: Longitude of the first point (C{degrees}).
341 @arg azi1: Azimuth at the first point (compass C{degrees}).
342 @arg s12: Distance between the points (C{meter}), can be negative.
343 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying
344 the quantities to be returned.
346 @return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2,
347 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1},
348 C{lon1}, C{azi1} and distance C{s12} always included.
350 @see: C++ U{GeodesicExact.Direct
351 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}
352 and Python U{Geodesic.Direct<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
353 '''
354 return self._GDictDirect(lat1, lon1, azi1, False, s12, outmask)
356 def Direct3(self, lat1, lon1, azi1, s12): # PYCHOK outmask
357 '''Return the destination lat, lon and reverse azimuth
358 (final bearing) in C{degrees}.
360 @return: L{Destination3Tuple}C{(lat, lon, final)}.
361 '''
362 r = self._GDictDirect(lat1, lon1, azi1, False, s12, Caps._AZIMUTH_LATITUDE_LONGITUDE)
363 return Destination3Tuple(r.lat2, r.lon2, r.azi2) # no iteration
365 def _DirectLine(self, ll1, azi12, s12=0, **caps_name):
366 '''(INTERNAL) Short-cut version.
367 '''
368 return self.DirectLine(ll1.lat, ll1.lon, azi12, s12, **caps_name)
370 def DirectLine(self, lat1, lon1, azi1, s12, caps=Caps.STANDARD, **name):
371 '''Define a L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as distance.
373 @arg lat1: Latitude of the first point (C{degrees}).
374 @arg lon1: Longitude of the first point (C{degrees}).
375 @arg azi1: Azimuth at the first point (compass C{degrees}).
376 @arg s12: Distance between the points (C{meter}), can be negative.
377 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying
378 the capabilities the L{GeodesicLineExact} instance
379 should possess, i.e., which quantities can be
380 returned by calls to L{GeodesicLineExact.Position}.
381 @kwarg name: Optional C{B{name}=NN} (C{str}).
383 @return: A L{GeodesicLineExact} instance.
385 @note: The third point of the L{GeodesicLineExact} is set to correspond
386 to the second point of the I{Inverse} geodesic problem.
388 @note: Latitude B{C{lat1}} should in the range C{[-90, +90]}.
390 @see: C++ U{GeodesicExact.DirectLine
391 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and
392 Python U{Geodesic.DirectLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
393 '''
394 return self._GenDirectLine(lat1, lon1, azi1, False, s12, caps, **name)
396 def _dn(self, sbet, cbet): # in gxline._GeodesicLineExact.__init__
397 '''(INTERNAL) Helper.
398 '''
399 if self.f < 0: # PYCHOK no cover
400 dn = sqrt(_1_0 - cbet**2 * self.e2) / self.f1
401 else:
402 dn = sqrt(_1_0 + sbet**2 * self.ep2)
403 return dn
405 @Property_RO
406 def e2(self):
407 '''Get the ellipsoid's I{(1st) eccentricity squared} (C{float}), M{f * (2 - f)}.
408 '''
409 return self.ellipsoid.e2
411 @Property_RO
412 def _e2a2(self):
413 '''(INTERNAL) Cache M{E.e2 * E.a2}.
414 '''
415 return self.e2 * self.ellipsoid.a2
417 @Property_RO
418 def _e2_f1(self):
419 '''(INTERNAL) Cache M{E.e2 * E.f1}.
420 '''
421 return self.e2 / self.f1
423 @Property_RO
424 def _eF(self):
425 '''(INTERNAL) Get the elliptic function, aka C{.E}.
426 '''
427 return _MODS.elliptic.Elliptic(k2=-self.ep2)
429 def _eF_reset_cHe2_f1(self, x, y):
430 '''(INTERNAL) Reset elliptic function and return M{cH * e2 / f1 * ...}.
431 '''
432 self._eF_reset_k2(x)
433 return y * self._eF.cH * self._e2_f1
435 def _eF_reset_k2(self, x):
436 '''(INTERNAL) Reset elliptic function and return C{k2}.
437 '''
438 ep2 = self.ep2
439 k2 = x**2 * ep2 # see .gxline._GeodesicLineExact._eF
440 self._eF.reset(k2=-k2, alpha2=-ep2) # kp2, alphap2 defaults
441 _update_glXs(self) # zap cached/memoized _GeodesicLineExact attrs
442 return k2
444 @Property_RO
445 def ellipsoid(self):
446 '''Get the ellipsoid (C{Ellipsoid}).
447 '''
448 return self.datum.ellipsoid
450 @Property_RO
451 def ep2(self):
452 '''Get the ellipsoid's I{2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)}.
453 '''
454 return self.ellipsoid.e22 # == self.e2 / self.f1**2
456 e22 = ep2 # for ellipsoid compatibility
458 @Property_RO
459 def _eTOL2(self):
460 '''(INTERNAL) The si12 threshold for "really short".
461 '''
462 return _eTOL2(self.f)
464 @Property_RO
465 def flattening(self):
466 '''Get the C{ellipsoid}'s I{flattening} (C{scalar}), M{(a - b) / a}, C{0} for spherical, negative for prolate.
467 '''
468 return self.ellipsoid.f
470 f = flattening
472 @Property_RO
473 def f1(self): # in .css.CassiniSoldner.reset
474 '''Get the C{ellipsoid}'s I{1 - flattening} (C{float}).
475 '''
476 return self.ellipsoid.f1
478 @Property_RO
479 def _f180(self):
480 '''(INTERNAL) Cached/memoized.
481 '''
482 return self.f * _180_0
484 def _GDictDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask=Caps.STANDARD):
485 '''(INTERNAL) As C{_GenDirect}, but returning a L{GDict}.
487 @return: A L{GDict} ...
488 '''
489 C = outmask if arcmode else (outmask | Caps.DISTANCE_IN)
490 glX = self.Line(lat1, lon1, azi1, C | Caps.LINE_OFF)
491 return glX._GDictPosition(arcmode, s12_a12, outmask)
493 def _GDictInverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD): # MCCABE 33, 41 vars
494 '''(INTERNAL) As C{_GenInverse}, but returning a L{GDict}.
496 @return: A L{GDict} ...
497 '''
498 Cs = Caps
499 if self._debug: # PYCHOK no cover
500 outmask |= Cs._DEBUG_INVERSE & self._debug
501 outmask &= Cs._OUT_MASK # incl. _SALPs_CALPs and _DEBUG_
502 # compute longitude difference carefully (with _diff182):
503 # result is in [-180, +180] but -180 is only for west-going
504 # geodesics, +180 is for east-going and meridional geodesics
505 lon12, lon12s = _diff182(lon1, lon2)
506 # see C{result} from geographiclib.geodesic.Inverse
507 if (outmask & Cs.LONG_UNROLL): # == (lon1 + lon12) + lon12s
508 r = GDict(lon1=lon1, lon2=fsumf_(lon1, lon12, lon12s))
509 else:
510 r = GDict(lon1=_norm180(lon1), lon2=_norm180(lon2))
511 if _K_2_0: # GeographicLib 2.0
512 # make longitude difference positive
513 lon12, lon_ = _unsigned2(lon12)
514 if lon_:
515 lon12s = -lon12s
516 lam12 = radians(lon12)
517 # calculate sincosd(_around(lon12 + correction))
518 slam12, clam12 = _sincos2de(lon12, lon12s)
519 # supplementary longitude difference
520 lon12s = fsumf_(_180_0, -lon12, -lon12s)
521 else: # GeographicLib 1.52
522 # make longitude difference positive and if very close
523 # to being on the same half-meridian, then make it so.
524 if lon12 < 0: # _signBit(lon12)
525 lon_, lon12 = True, -_around(lon12)
526 lon12s = _around(fsumf_(_180_0, -lon12, lon12s))
527 else:
528 lon_, lon12 = False, _around(lon12)
529 lon12s = _around(fsumf_(_180_0, -lon12, -lon12s))
530 lam12 = radians(lon12)
531 if lon12 > _90_0:
532 slam12, clam12 = _sincos2d(lon12s)
533 clam12 = -clam12
534 else:
535 slam12, clam12 = _sincos2(lam12)
536 # If really close to the equator, treat as on equator.
537 lat1 = _around(_fix90(lat1))
538 lat2 = _around(_fix90(lat2))
539 r.set_(lat1=lat1, lat2=lat2)
540 # Swap points so that point with higher (abs) latitude is
541 # point 1. If one latitude is a NAN, then it becomes lat1.
542 swap_ = fabs(lat1) < fabs(lat2) or isnan(lat2)
543 if swap_:
544 lat1, lat2 = lat2, lat1
545 lon_ = not lon_
546 if _signBit(lat1):
547 lat_ = False # note, False
548 else: # make lat1 <= -0
549 lat_ = True # note, True
550 lat1, lat2 = -lat1, -lat2
551 # Now 0 <= lon12 <= 180, -90 <= lat1 <= -0 and lat1 <= lat2 <= -lat1
552 # and lat_, lon_, swap_ register the transformation to bring the
553 # coordinates to this canonical form, where False means no change
554 # made. We make these transformations so that there are few cases
555 # to check, e.g., on verifying quadrants in atan2. In addition,
556 # this enforces some symmetries in the results returned.
558 # Initialize for the meridian. No longitude calculation is
559 # done in this case to let the parameter default to 0.
560 sbet1, cbet1 = self._sinf1cos2d(lat1)
561 sbet2, cbet2 = self._sinf1cos2d(lat2)
562 # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure
563 # of the |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1),
564 # abs(sbet2) + sbet1 is a better measure. This logic is used
565 # in assigning calp2 in _Lambda6. Sometimes these quantities
566 # vanish and in that case we force bet2 = +/- bet1 exactly. An
567 # example where is is necessary is the inverse problem
568 # 48.522876735459 0 -48.52287673545898293 179.599720456223079643
569 # which failed with Visual Studio 10 (Release and Debug)
570 if cbet1 < -sbet1:
571 if cbet2 == cbet1:
572 sbet2 = copysign(sbet1, sbet2)
573 elif fabs(sbet2) == -sbet1:
574 cbet2 = cbet1
576 p = _PDict(sbet1=sbet1, cbet1=cbet1, dn1=self._dn(sbet1, cbet1),
577 sbet2=sbet2, cbet2=cbet2, dn2=self._dn(sbet2, cbet2))
579 _meridian = _b = True # i.e. not meridian, not b
580 if lat1 == -90 or slam12 == 0:
581 # Endpoints are on a single full meridian,
582 # so the geodesic might lie on a meridian.
583 salp1, calp1 = slam12, clam12 # head to target lon
584 salp2, calp2 = _0_0, _1_0 # then head north
585 # tan(bet) = tan(sig) * cos(alp)
586 p.set_sigs(sbet1, calp1 * cbet1, sbet2, calp2 * cbet2)
587 # sig12 = sig2 - sig1
588 sig12 = _atan12(sbet1, p.csig1, sbet2, p.csig2, sineg0=True) # PYCHOK csig*
589 s12x, m12x, _, \
590 M12, M21 = self._Length5(sig12, outmask | Cs.REDUCEDLENGTH, p)
591 # Add the check for sig12 since zero length geodesics
592 # might yield m12 < 0. Test case was
593 # echo 20.001 0 20.001 0 | GeodSolve -i
594 # In fact, we will have sig12 > PI/2 for meridional
595 # geodesic which is not a shortest path.
596 if m12x >= 0 or sig12 < _1_0:
597 # Need at least 2 to handle 90 0 90 180
598 # Prevent negative s12 or m12 from geographiclib 1.52
599 if sig12 < _TINY3 or (sig12 < _TOL0 and (s12x < 0 or m12x < 0)):
600 sig12 = m12x = s12x = _0_0
601 else:
602 _b = False # apply .b to s12x, m12x
603 _meridian = False
604 C = 1
605 # else: # m12 < 0, prolate and too close to anti-podal
606 # _meridian = True
607 a12 = _0_0 # if _b else degrees(sig12)
609 if _meridian:
610 _b = sbet1 == 0 and (self.f <= 0 or lon12s >= self._f180) # and sbet2 == 0
611 if _b: # Geodesic runs along equator
612 calp1 = calp2 = _0_0
613 salp1 = salp2 = _1_0
614 sig12 = lam12 / self.f1 # == omg12
615 somg12, comg12 = _sincos2(sig12)
616 m12x = self.b * somg12
617 s12x = self.a * lam12
618 M12 = M21 = comg12
619 a12 = lon12 / self.f1
620 C = 2
621 else:
622 # Now point1 and point2 belong within a hemisphere bounded by a
623 # meridian and geodesic is neither meridional or equatorial.
624 p.set_(slam12=slam12, clam12=clam12)
625 # Figure a starting point for Newton's method
626 sig12, salp1, calp1, \
627 salp2, calp2, dnm = self._InverseStart6(lam12, p)
628 if sig12 is None: # use Newton's method
629 # pre-compute the constant _Lambda6 term, once
630 p.set_(bet12=None if cbet2 == cbet1 and fabs(sbet2) == -sbet1 else
631 (((cbet1 + cbet2) * (cbet2 - cbet1)) if cbet1 < -sbet1 else
632 ((sbet1 + sbet2) * (sbet1 - sbet2))))
633 sig12, salp1, calp1, \
634 salp2, calp2, domg12 = self._Newton6(salp1, calp1, p)
635 s12x, m12x, _, M12, M21 = self._Length5(sig12, outmask, p)
636 if (outmask & Cs.AREA):
637 # omg12 = lam12 - domg12
638 s, c = _sincos2(domg12)
639 somg12, comg12 = _sincos12(s, c, slam12, clam12)
640 C = 3 # Newton
641 else: # from _InverseStart6: dnm, salp*, calp*
642 C = 4 # Short lines
643 s, c = _sincos2(sig12 / dnm)
644 m12x = dnm**2 * s
645 s12x = dnm * sig12
646 M12 = M21 = c
647 if (outmask & Cs.AREA):
648 somg12, comg12 = _sincos2(lam12 / (self.f1 * dnm))
650 else: # _meridian is False
651 somg12 = comg12 = NAN
653 r.set_(a12=a12 if _b else degrees(sig12)) # in [0, 180]
655 if (outmask & Cs.DISTANCE):
656 r.set_(s12=unsigned0(s12x if _b else (self.b * s12x)))
658 if (outmask & Cs.REDUCEDLENGTH):
659 r.set_(m12=unsigned0(m12x if _b else (self.b * m12x)))
661 if (outmask & Cs.GEODESICSCALE):
662 if swap_:
663 M12, M21 = M21, M12
664 r.set_(M12=unsigned0(M12),
665 M21=unsigned0(M21))
667 if (outmask & Cs.AREA):
668 S12 = self._InverseArea(_meridian, salp1, calp1,
669 salp2, calp2,
670 somg12, comg12, p)
671 if _xor(swap_, lat_, lon_):
672 S12 = -S12
673 r.set_(S12=unsigned0(S12))
675 if (outmask & (Cs.AZIMUTH | Cs._SALPs_CALPs)):
676 if swap_:
677 salp1, salp2 = salp2, salp1
678 calp1, calp2 = calp2, calp1
679 if _xor(swap_, lon_):
680 salp1, salp2 = -salp1, -salp2
681 if _xor(swap_, lat_):
682 calp1, calp2 = -calp1, -calp2
684 if (outmask & Cs.AZIMUTH):
685 r.set_(azi1=_atan2d(salp1, calp1),
686 azi2=_atan2d_reverse(salp2, calp2, reverse=outmask & Cs.REVERSE2))
687 if (outmask & Cs._SALPs_CALPs):
688 r.set_(salp1=salp1, calp1=calp1,
689 salp2=salp2, calp2=calp2)
691 if (outmask & Cs._DEBUG_INVERSE): # PYCHOK no cover
692 E, eF = self.ellipsoid, self._eF
693 p.set_(C=C, a=self.a, f=self.f, f1=self.f1,
694 e=E.e, e2=self.e2, ep2=self.ep2,
695 c2=E.c2, c2x=self.c2x,
696 eFcD=eF.cD, eFcE=eF.cE, eFcH=eF.cH,
697 eFk2=eF.k2, eFa2=eF.alpha2)
698 p.update(r) # r overrides p
699 r = p.toGDict()
700 return self._iter2tion(r, **p)
702 def _GenDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask=Caps.STANDARD):
703 '''(INTERNAL) The general I{Inverse} geodesic calculation.
705 @return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2,
706 s12, m12, M12, M21, S12)}.
707 '''
708 r = self._GDictDirect(lat1, lon1, azi1, arcmode, s12_a12, outmask)
709 return r.toDirect9Tuple()
711 def _GenDirectLine(self, lat1, lon1, azi1, arcmode, s12_a12, caps, **name):
712 '''(INTERNAL) Helper for C{ArcDirectLine} and C{DirectLine}.
714 @return: A L{GeodesicLineExact} instance.
715 '''
716 azi1 = _norm180(azi1)
717 # guard against underflow in salp0. Also -0 is converted to +0.
718 s, c = _sincos2d(_around(azi1))
719 C = caps if arcmode else (caps | Caps.DISTANCE_IN)
720 return _GeodesicLineExact(self, lat1, lon1, azi1, C,
721 self._debug, s, c, **name)._GenSet(arcmode, s12_a12)
723 def _GenInverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD):
724 '''(INTERNAL) The general I{Inverse} geodesic calculation.
726 @return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2,
727 m12, M12, M21, S12)}.
728 '''
729 r = self._GDictInverse(lat1, lon1, lat2, lon2, outmask | Caps._SALPs_CALPs)
730 return r.toInverse10Tuple()
732 def _Inverse(self, ll1, ll2, wrap, **outmask):
733 '''(INTERNAL) Short-cut version, see .base.ellipsoidalDI.intersecant2.
734 '''
735 if wrap:
736 ll2 = _unrollon(ll1, _Wrap.point(ll2))
737 return self.Inverse(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **outmask)
739 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD):
740 '''Perform the I{Inverse} geodesic calculation.
742 @arg lat1: Latitude of the first point (C{degrees}).
743 @arg lon1: Longitude of the first point (C{degrees}).
744 @arg lat2: Latitude of the second point (C{degrees}).
745 @arg lon2: Longitude of the second point (C{degrees}).
746 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying
747 the quantities to be returned.
749 @return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2,
750 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1},
751 C{lon1}, C{azi1} and distance C{s12} always included.
753 @note: The third point of the L{GeodesicLineExact} is set to correspond
754 to the second point of the I{Inverse} geodesic problem.
756 @note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}.
758 @see: C++ U{GeodesicExact.InverseLine
759 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and
760 Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
761 '''
762 return self._GDictInverse(lat1, lon1, lat2, lon2, outmask)
764 def Inverse1(self, lat1, lon1, lat2, lon2, wrap=False):
765 '''Return the non-negative, I{angular} distance in C{degrees}.
767 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
768 B{C{lat2}} and B{C{lon2}} (C{bool}).
769 '''
770 # see .FrechetKarney.distance, .HausdorffKarney._distance
771 # and .HeightIDWkarney._distances
772 if wrap:
773 _, lat2, lon2 = _Wrap.latlon3(lat1, lat2, lon2, True) # _Geodesic.LONG_UNROLL
774 return fabs(self._GDictInverse(lat1, lon1, lat2, lon2, Caps._ANGLE_ONLY).a12)
776 def Inverse3(self, lat1, lon1, lat2, lon2): # PYCHOK outmask
777 '''Return the distance in C{meter} and the forward and
778 reverse azimuths (initial and final bearing) in C{degrees}.
780 @return: L{Distance3Tuple}C{(distance, initial, final)}.
781 '''
782 r = self._GDictInverse(lat1, lon1, lat2, lon2, Caps.AZIMUTH_DISTANCE)
783 return Distance3Tuple(r.s12, wrap360(r.azi1), wrap360(r.azi2),
784 iteration=r.iteration)
786 def _InverseLine(self, ll1, ll2, wrap, **caps_name):
787 '''(INTERNAL) Short-cut version.
788 '''
789 if wrap:
790 ll2 = _unrollon(ll1, _Wrap.point(ll2))
791 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **caps_name)
793 def InverseLine(self, lat1, lon1, lat2, lon2, caps=Caps.STANDARD, **name):
794 '''Define a L{GeodesicLineExact} in terms of the I{Inverse} geodesic problem.
796 @arg lat1: Latitude of the first point (C{degrees}).
797 @arg lon1: Longitude of the first point (C{degrees}).
798 @arg lat2: Latitude of the second point (C{degrees}).
799 @arg lon2: Longitude of the second point (C{degrees}).
800 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying
801 the capabilities the L{GeodesicLineExact} instance
802 should possess, i.e., which quantities can be
803 returned by calls to L{GeodesicLineExact.Position}
804 and L{GeodesicLineExact.ArcPosition}.
805 @kwarg name: Optional C{B{name}=NN} (C{str}).
807 @return: A L{GeodesicLineExact} instance.
809 @note: The third point of the L{GeodesicLineExact} is set to correspond
810 to the second point of the I{Inverse} geodesic problem.
812 @note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}.
814 @see: C++ U{GeodesicExact.InverseLine
815 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and
816 Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
817 '''
818 Cs = Caps
819 r = self._GDictInverse(lat1, lon1, lat2, lon2, Cs._SALPs_CALPs) # No need for AZIMUTH
820 C = (caps | Cs.DISTANCE) if (caps & (Cs.DISTANCE_IN & Cs._OUT_MASK)) else caps
821 azi1 = _atan2d(r.salp1, r.calp1)
822 return _GeodesicLineExact(self, lat1, lon1, azi1, C, # ensure a12 is distance
823 self._debug, r.salp1, r.calp1, **name)._GenSet(True, r.a12)
825 def _InverseArea(self, _meridian, salp1, calp1, # PYCHOK 9 args
826 salp2, calp2,
827 somg12, comg12, p):
828 '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length.
830 @return: Area C{S12}.
831 '''
832 # from _Lambda6: sin(alp1) * cos(bet1) = sin(alp0), calp0 > 0
833 salp0, calp0 = _sin1cos2(salp1, calp1, p.sbet1, p.cbet1)
834 A4 = calp0 * salp0
835 if A4:
836 # from _Lambda6: tan(bet) = tan(sig) * cos(alp)
837 k2 = calp0**2 * self.ep2
838 C4a = self._C4f_k2(k2)
839 B41 = _cosSeries(C4a, *_norm2(p.sbet1, calp1 * p.cbet1))
840 B42 = _cosSeries(C4a, *_norm2(p.sbet2, calp2 * p.cbet2))
841 # multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
842 A4 *= self._e2a2
843 S12 = A4 * (B42 - B41)
844 else: # avoid problems with indeterminate sig1, sig2 on equator
845 A4 = B41 = B42 = k2 = S12 = _0_0
847 if (_meridian and # omg12 < 3/4 * PI
848 comg12 > -_SQRT2_2 and # lon diff not too big
849 (p.sbet2 - p.sbet1) < _1_75): # lat diff not too big
850 # use tan(Gamma/2) = tan(omg12/2) *
851 # (tan(bet1/2) + tan(bet2/2)) /
852 # (tan(bet1/2) * tan(bet2/2) + 1))
853 # with tan(x/2) = sin(x) / (1 + cos(x))
854 dbet1 = p.cbet1 + _1_0
855 dbet2 = p.cbet2 + _1_0
856 domg12 = comg12 + _1_0
857 salp12 = (p.sbet1 * dbet2 + dbet1 * p.sbet2) * somg12
858 calp12 = (p.sbet1 * p.sbet2 + dbet1 * dbet2) * domg12
859 alp12 = atan2(salp12, calp12) * _2_0
860 else:
861 # alp12 = alp2 - alp1, used in atan2, no need to normalize
862 salp12, calp12 = _sincos12(salp1, calp1, salp2, calp2)
863 # The right thing appears to happen if alp1 = +/-180 and
864 # alp2 = 0, viz salp12 = -0 and alp12 = -180. However,
865 # this depends on the sign being attached to 0 correctly.
866 # Following ensures the correct behavior.
867 if salp12 == 0 and calp12 < 0:
868 alp12 = _copysign(PI, calp1)
869 else:
870 alp12 = atan2(salp12, calp12)
872 p.set_(alp12=alp12, A4=A4, B41=B41, B42=B42, k2=k2)
873 return S12 + self.c2x * alp12
875 def _InverseStart6(self, lam12, p):
876 '''(INTERNAL) Return a starting point for Newton's method in
877 C{salp1} and C{calp1} indicated by C{sig12=None}. If
878 Newton's method doesn't need to be used, return also
879 C{salp2}, C{calp2}, C{dnm} and C{sig12} non-C{None}.
881 @return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, dnm)}
882 and C{p.set_sigs} updated for Newton, C{sig12=None}.
883 '''
884 sig12 = None # use Newton
885 salp1 = calp1 = salp2 = calp2 = dnm = NAN
887 # bet12 = bet2 - bet1 in [0, PI)
888 sbet12, cbet12 = _sincos12(p.sbet1, p.cbet1, p.sbet2, p.cbet2)
889 shortline = cbet12 >= 0 and sbet12 < _0_5 and (p.cbet2 * lam12) < _0_5
890 if shortline:
891 # sin((bet1 + bet2)/2)^2 = (sbet1 + sbet2)^2 / (
892 # (cbet1 + cbet2)^2 + (sbet1 + sbet2)^2)
893 t = (p.sbet1 + p.sbet2)**2
894 s = t / ((p.cbet1 + p.cbet2)**2 + t)
895 dnm = sqrt(_1_0 + self.ep2 * s)
896 somg12, comg12 = _sincos2(lam12 / (self.f1 * dnm))
897 else:
898 somg12, comg12 = p.slam12, p.clam12
900 # bet12a = bet2 + bet1 in (-PI, 0], note -sbet1
901 sbet12a, cbet12a = _sincos12(-p.sbet1, p.cbet1, p.sbet2, p.cbet2)
903 c = fabs(comg12) + _1_0 # == (1 - comg12) if comg12 < 0
904 s = somg12**2 / c
905 t = p.sbet1 * p.cbet2 * s
906 salp1 = p.cbet2 * somg12
907 calp1 = (sbet12a - t) if comg12 < 0 else (sbet12 + t)
909 ssig12 = _hypot(salp1, calp1)
910 csig12 = p.sbet1 * p.sbet2 + p.cbet1 * p.cbet2 * comg12
912 if shortline and ssig12 < self._eTOL2: # really short lines
913 t = c if comg12 < 0 else s
914 salp2, calp2 = _norm2(somg12 * p.cbet1,
915 sbet12 - p.cbet1 * p.sbet2 * t)
916 sig12 = atan2(ssig12, csig12) # do not use Newton
918 elif (self._n_0_1 or # Skip astroid calc if too eccentric
919 csig12 >= 0 or ssig12 >= (p.cbet1**2 * self._n6PI)):
920 pass # nothing to do, 0th order spherical approximation OK
922 else:
923 # Scale lam12 and bet2 to x, y coordinate system where antipodal
924 # point is at origin and singular point is at y = 0, x = -1
925 lam12x = atan2(-p.slam12, -p.clam12) # lam12 - PI
926 f = self.f
927 if f < 0: # PYCHOK no cover
928 # ssig1=sbet1, csig1=-cbet1, ssig2=sbet2, csig2=cbet2
929 p.set_sigs(p.sbet1, -p.cbet1, p.sbet2, p.cbet2)
930 # if lon12 = 180, this repeats a calculation made in Inverse
931 _, m12b, m0, _, _ = self._Length5(atan2(sbet12a, cbet12a) + PI,
932 Caps.REDUCEDLENGTH, p)
933 t = p.cbet1 * PI # x = dlat, y = dlon
934 x = m12b / (t * p.cbet2 * m0) - _1_0
935 sca = (sbet12a / (x * p.cbet1)) if x < -_0_01 else (-f * t)
936 y = lam12x / sca
937 else: # f >= 0, however f == 0 does not get here
938 sca = self._eF_reset_cHe2_f1(p.sbet1, p.cbet1 * _2_0)
939 x = lam12x / sca # dlon
940 y = sbet12a / (sca * p.cbet1) # dlat
942 if y > _TOL1 and x > -_THR1: # strip near cut
943 if f < 0: # PYCHOK no cover
944 calp1 = max( _0_0, x) if x > _TOL1 else max(_N_1_0, x)
945 salp1 = sqrt(_1_0 - calp1**2)
946 else:
947 salp1 = min( _1_0, -x)
948 calp1 = -sqrt(_1_0 - salp1**2)
949 else:
950 # Estimate alp1, by solving the astroid problem.
951 #
952 # Could estimate alpha1 = theta + PI/2, directly, i.e.,
953 # calp1 = y/k; salp1 = -x/(1+k); for _f >= 0
954 # calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check)
955 #
956 # However, it's better to estimate omg12 from astroid and use
957 # spherical formula to compute alp1. This reduces the mean
958 # number of Newton iterations for astroid cases from 2.24
959 # (min 0, max 6) to 2.12 (min 0, max 5). The changes in the
960 # number of iterations are as follows:
961 #
962 # change percent
963 # 1 5
964 # 0 78
965 # -1 16
966 # -2 0.6
967 # -3 0.04
968 # -4 0.002
969 #
970 # The histogram of iterations is (m = number of iterations
971 # estimating alp1 directly, n = number of iterations
972 # estimating via omg12, total number of trials = 148605):
973 #
974 # iter m n
975 # 0 148 186
976 # 1 13046 13845
977 # 2 93315 102225
978 # 3 36189 32341
979 # 4 5396 7
980 # 5 455 1
981 # 6 56 0
982 #
983 # omg12 is near PI, estimate work with omg12a = PI - omg12
984 k = _Astroid(x, y)
985 sca *= (y * (k + _1_0) / k) if f < 0 else \
986 (x * k / (k + _1_0))
987 s, c = _sincos2(-sca) # omg12a
988 # update spherical estimate of alp1 using omg12 instead of lam12
989 salp1 = p.cbet2 * s
990 calp1 = sbet12a - s * salp1 * p.sbet1 / (c + _1_0) # c = -c
992 # sanity check on starting guess. Backwards check allows NaN through.
993 salp1, calp1 = _norm2(salp1, calp1) if salp1 > 0 else (_1_0, _0_0)
995 return sig12, salp1, calp1, salp2, calp2, dnm
997 def _Lambda6(self, salp1, calp1, diffp, p):
998 '''(INTERNAL) Helper.
1000 @return: 6-Tuple C{(lam12, sig12, salp2, calp2, domg12,
1001 dlam12} and C{p.set_sigs} updated.
1002 '''
1003 eF = self._eF
1004 f1 = self.f1
1006 if p.sbet1 == calp1 == 0: # PYCHOK no cover
1007 # Break degeneracy of equatorial line
1008 calp1 = -_TINY
1010 # sin(alp1) * cos(bet1) = sin(alp0), # calp0 > 0
1011 salp0, calp0 = _sin1cos2(salp1, calp1, p.sbet1, p.cbet1)
1012 # tan(bet1) = tan(sig1) * cos(alp1)
1013 # tan(omg1) = sin(alp0) * tan(sig1)
1014 # = sin(bet1) * tan(alp1)
1015 somg1 = salp0 * p.sbet1
1016 comg1 = calp1 * p.cbet1
1017 ssig1, csig1 = _norm2(p.sbet1, comg1)
1018 # Without normalization we have schi1 = somg1
1019 cchi1 = f1 * p.dn1 * comg1
1021 # Enforce symmetries in the case abs(bet2) = -bet1.
1022 # Need to be careful about this case, since this can
1023 # yield singularities in the Newton iteration.
1024 # sin(alp2) * cos(bet2) = sin(alp0)
1025 salp2 = (salp0 / p.cbet2) if p.cbet2 != p.cbet1 else salp1
1026 # calp2 = sqrt(1 - sq(salp2))
1027 # = sqrt(sq(calp0) - sq(sbet2)) / cbet2
1028 # and subst for calp0 and rearrange to give (choose
1029 # positive sqrt to give alp2 in [0, PI/2]).
1030 calp2 = fabs(calp1) if p.bet12 is None else (
1031 sqrt((calp1 * p.cbet1)**2 + p.bet12) / p.cbet2)
1032 # tan(bet2) = tan(sig2) * cos(alp2)
1033 # tan(omg2) = sin(alp0) * tan(sig2).
1034 somg2 = salp0 * p.sbet2
1035 comg2 = calp2 * p.cbet2
1036 ssig2, csig2 = _norm2(p.sbet2, comg2)
1037 # without normalization we have schi2 = somg2
1038 cchi2 = f1 * p.dn2 * comg2
1040 # omg12 = omg2 - omg1, limit to [0, PI]
1041 somg12, comg12 = _sincos12(somg1, comg1, somg2, comg2, sineg0=True)
1042 # chi12 = chi2 - chi1, limit to [0, PI]
1043 schi12, cchi12 = _sincos12(somg1, cchi1, somg2, cchi2, sineg0=True)
1045 p.set_sigs(ssig1, csig1, ssig2, csig2)
1046 # sig12 = sig2 - sig1, limit to [0, PI]
1047 sig12 = _atan12(ssig1, csig1, ssig2, csig2, sineg0=True)
1049 eta12 = self._eF_reset_cHe2_f1(calp0, salp0) * _2__PI # then ...
1050 eta12 *= fsum1f_(eF.deltaH(*p.sncndn2),
1051 -eF.deltaH(*p.sncndn1), sig12)
1052 # eta = chi12 - lam12
1053 lam12 = _atan12(p.slam12, p.clam12, schi12, cchi12) - eta12
1054 # domg12 = chi12 - omg12 - deta12
1055 domg12 = _atan12( somg12, comg12, schi12, cchi12) - eta12
1057 dlam12 = NAN # dv > 0 in ._Newton6
1058 if diffp:
1059 d = calp2 * p.cbet2
1060 if d:
1061 _, dlam12, _, _, _ = self._Length5(sig12, Caps.REDUCEDLENGTH, p)
1062 dlam12 *= f1 / d
1063 elif p.sbet1:
1064 dlam12 = -f1 * p.dn1 * _2_0 / p.sbet1
1066 # p.set_(deta12=-eta12, lam12=lam12)
1067 return lam12, sig12, salp2, calp2, domg12, dlam12
1069 def _Length5(self, sig12, outmask, p):
1070 '''(INTERNAL) Return M{m12b = (reduced length) / self.b} and
1071 calculate M{s12b = distance / self.b} and M{m0}, the
1072 coefficient of secular term in expression for reduced
1073 length and the geodesic scales C{M12} and C{M21}.
1075 @return: 5-Tuple C{(s12b, m12b, m0, M12, M21)}.
1076 '''
1077 s12b = m12b = m0 = M12 = M21 = NAN
1079 Cs = Caps
1080 eF = self._eF
1082 # outmask &= Cs._OUT_MASK
1083 if (outmask & Cs.DISTANCE):
1084 # Missing a factor of self.b
1085 s12b = eF.cE * _2__PI * fsum1f_(eF.deltaE(*p.sncndn2),
1086 -eF.deltaE(*p.sncndn1), sig12)
1088 if (outmask & Cs._REDUCEDLENGTH_GEODESICSCALE):
1089 m0x = -eF.k2 * eF.cD * _2__PI
1090 J12 = -m0x * fsum1f_(eF.deltaD(*p.sncndn2),
1091 -eF.deltaD(*p.sncndn1), sig12)
1092 if (outmask & Cs.REDUCEDLENGTH):
1093 m0 = m0x
1094 # Missing a factor of self.b. Add parens around
1095 # (csig1 * ssig2) and (ssig1 * csig2) to ensure
1096 # accurate cancellation for coincident points.
1097 m12b = fsum1f_(p.dn2 * (p.csig1 * p.ssig2),
1098 -p.dn1 * (p.ssig1 * p.csig2),
1099 J12 * (p.csig1 * p.csig2))
1100 if (outmask & Cs.GEODESICSCALE):
1101 M12 = M21 = p.ssig1 * p.ssig2 + \
1102 p.csig1 * p.csig2
1103 t = (p.cbet1 - p.cbet2) * self.ep2 * \
1104 (p.cbet1 + p.cbet2) / (p.dn1 + p.dn2)
1105 M12 += (p.ssig2 * t + p.csig2 * J12) * p.ssig1 / p.dn1
1106 M21 -= (p.ssig1 * t + p.csig1 * J12) * p.ssig2 / p.dn2
1108 return s12b, m12b, m0, M12, M21
1110 def Line(self, lat1, lon1, azi1, caps=Caps.ALL, **name):
1111 '''Set up a L{GeodesicLineExact} to compute several points
1112 on a single geodesic.
1114 @arg lat1: Latitude of the first point (C{degrees}).
1115 @arg lon1: Longitude of the first point (C{degrees}).
1116 @arg azi1: Azimuth at the first point (compass C{degrees}).
1117 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying
1118 the capabilities the L{GeodesicLineExact} instance
1119 should possess, i.e., which quantities can be
1120 returnedby calls to L{GeodesicLineExact.Position}
1121 and L{GeodesicLineExact.ArcPosition}.
1122 @kwarg name: Optional C{B{name}=NN} (C{str}).
1124 @return: A L{GeodesicLineExact} instance.
1126 @note: If the point is at a pole, the azimuth is defined by keeping
1127 B{C{lon1}} fixed, writing C{B{lat1} = ±(90 − ε)}, and taking
1128 the limit C{ε → 0+}.
1130 @see: C++ U{GeodesicExact.Line
1131 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}
1132 and Python U{Geodesic.Line<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.
1133 '''
1134 return _GeodesicLineExact(self, lat1, lon1, azi1, caps, self._debug, **name)
1136 @Property_RO
1137 def n(self):
1138 '''Get the C{ellipsoid}'s I{3rd flattening} (C{scalar}), M{f / (2 - f) == (a - b) / (a + b)}.
1139 '''
1140 return self.ellipsoid.n
1142 @Property_RO
1143 def _n_0_1(self):
1144 '''(INTERNAL) Cached once.
1145 '''
1146 return fabs(self.n) > _0_1
1148 @Property_RO
1149 def _n6PI(self):
1150 '''(INTERNAL) Cached once.
1151 '''
1152 return fabs(self.n) * _6_0 * PI
1154 def _Newton6(self, salp1, calp1, p):
1155 '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length.
1157 @return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, domg12)}
1158 and C{p.iter} and C{p.trip} updated.
1159 '''
1160 _abs = fabs
1161 # This is a straightforward solution of f(alp1) = lambda12(alp1) -
1162 # lam12 = 0 with one wrinkle. f(alp) has exactly one root in the
1163 # interval (0, PI) and its derivative is positive at the root.
1164 # Thus f(alp) is positive for alp > alp1 and negative for alp < alp1.
1165 # During the course of the iteration, a range (alp1a, alp1b) is
1166 # maintained which brackets the root and with each evaluation of
1167 # f(alp) the range is shrunk, if possible. Newton's method is
1168 # restarted whenever the derivative of f is negative (because the
1169 # new value of alp1 is then further from the solution) or if the
1170 # new estimate of alp1 lies outside (0,PI); in this case, the new
1171 # starting guess is taken to be (alp1a + alp1b) / 2.
1172 salp1a = salp1b = _TINY
1173 calp1a, calp1b = _1_0, _N_1_0
1174 MAXIT1, TOL0 = _MAXIT1, _TOL0
1175 HALF, TOLb = _0_5, _TOLb
1176 tripb, TOLv = False, TOL0
1177 for i in range(_MAXIT2):
1178 # 1/4 meridian = 10e6 meter and random input,
1179 # estimated max error in nm (nano meter, by
1180 # checking Inverse problem by Direct).
1181 #
1182 # max iterations
1183 # log2(b/a) error mean sd
1184 # -7 387 5.33 3.68
1185 # -6 345 5.19 3.43
1186 # -5 269 5.00 3.05
1187 # -4 210 4.76 2.44
1188 # -3 115 4.55 1.87
1189 # -2 69 4.35 1.38
1190 # -1 36 4.05 1.03
1191 # 0 15 0.01 0.13
1192 # 1 25 5.10 1.53
1193 # 2 96 5.61 2.09
1194 # 3 318 6.02 2.74
1195 # 4 985 6.24 3.22
1196 # 5 2352 6.32 3.44
1197 # 6 6008 6.30 3.45
1198 # 7 19024 6.19 3.30
1199 v, sig12, salp2, calp2, \
1200 domg12, dv = self._Lambda6(salp1, calp1, i < MAXIT1, p)
1202 # 2 * _TOL0 is approximately 1 ulp [0, PI]
1203 # reversed test to allow escape with NaNs
1204 if tripb or _abs(v) < TOLv:
1205 break
1206 # update bracketing values
1207 if v > 0 and (i > MAXIT1 or (calp1 / salp1) > (calp1b / salp1b)):
1208 salp1b, calp1b = salp1, calp1
1209 elif v < 0 and (i > MAXIT1 or (calp1 / salp1) < (calp1a / salp1a)):
1210 salp1a, calp1a = salp1, calp1
1212 if i < MAXIT1 and dv > 0:
1213 dalp1 = -v / dv
1214 if _abs(dalp1) < PI:
1215 s, c = _sincos2(dalp1)
1216 # nalp1 = alp1 + dalp1
1217 s, c = _sincos12(-s, c, salp1, calp1)
1218 if s > 0:
1219 salp1, calp1 = _norm2(s, c)
1220 # in some regimes we don't get quadratic convergence
1221 # because slope -> 0. So use convergence conditions
1222 # based on epsilon instead of sqrt(epsilon)
1223 TOLv = TOL0 if _abs(v) > _TOL016 else _TOL08
1224 continue
1225 TOLv = TOL0
1226 # Either dv was not positive or updated value was outside
1227 # legal range. Use the midpoint of the bracket as the next
1228 # estimate. This mechanism is not needed for the WGS84
1229 # ellipsoid, but it does catch problems with more eccentric
1230 # ellipsoids. Its efficacy is such for the WGS84 test set
1231 # with the starting guess set to alp1 = 90 deg: the WGS84
1232 # test set: mean = 5.21, stdev = 3.93, max = 24 and WGS84
1233 # with random input: mean = 4.74, stdev = 0.99
1234 salp1, calp1 = _norm2((salp1a + salp1b) * HALF,
1235 (calp1a + calp1b) * HALF)
1236 tripb = fsum1f_(calp1a, -calp1, _abs(salp1a - salp1)) < TOLb or \
1237 fsum1f_(calp1b, -calp1, _abs(salp1b - salp1)) < TOLb
1238 else:
1239 raise GeodesicError(Fmt.no_convergence(v, TOLv), txt=repr(self)) # self.toRepr()
1241 p.set_(iter=i, trip=tripb) # like .geodsolve._GDictInvoke: iter NOT iteration!
1242 return sig12, salp1, calp1, salp2, calp2, domg12
1244 Polygon = Area # for C{geographiclib} compatibility
1246 def _sinf1cos2d(self, lat):
1247 '''(INTERNAL) Helper, also for C{_G_GeodesicLineExact}.
1248 '''
1249 sbet, cbet = _sincos2d(lat)
1250 # ensure cbet1 = +epsilon at poles; doing the fix on beta means
1251 # that sig12 will be <= 2*tiny for two points at the same pole
1252 sbet, cbet = _norm2(sbet * self.f1, cbet)
1253 return sbet, (cbet if fabs(cbet) > _TINY else _TINY)
1255 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature
1256 '''Return this C{GeodesicExact} as string.
1258 @kwarg prec: The C{float} precision, number of decimal digits (0..9).
1259 Trailing zero decimals are stripped for B{C{prec}} values
1260 of 1 and above, but kept for negative B{C{prec}} values.
1261 @kwarg sep: Separator to join (C{str}).
1263 @return: Tuple items (C{str}).
1264 '''
1265 d = dict(ellipsoid=self.ellipsoid, C4order=self.C4order)
1266 return sep.join(pairs(d, prec=prec))
1269class GeodesicLineExact(_GeodesicLineExact):
1270 '''A pure Python version of I{Karney}'s C++ class U{GeodesicLineExact
1271 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicLineExact.html>},
1272 modeled after I{Karney}'s Python class U{geodesicline.GeodesicLine<https://GitHub.com/
1273 geographiclib/geographiclib-python>}.
1274 '''
1276 def __init__(self, geodesic, lat1, lon1, azi1, caps=Caps.STANDARD, **name):
1277 '''New L{GeodesicLineExact} instance, allowing points to be found along
1278 a geodesic starting at C{(B{lat1}, B{lon1})} with azimuth B{C{azi1}}.
1280 @arg geodesic: The geodesic to use (L{GeodesicExact}).
1281 @arg lat1: Latitude of the first point (C{degrees}).
1282 @arg lon1: Longitude of the first point (C{degrees}).
1283 @arg azi1: Azimuth at the first points (compass C{degrees}).
1284 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying
1285 the capabilities the L{GeodesicLineExact} instance
1286 should possess, i.e., which quantities can be
1287 returned by calls to L{GeodesicLineExact.Position}
1288 and L{GeodesicLineExact.ArcPosition}.
1289 @kwarg name: Optional C{B{name}=NN} (C{str}).
1291 @raise TypeError: Invalid B{C{geodesic}}.
1292 '''
1293 _xinstanceof(GeodesicExact, geodesic=geodesic)
1294 if (caps & Caps.LINE_OFF): # copy to avoid updates
1295 geodesic = geodesic.copy(deep=False, name=_UNDER_(NN, geodesic.name)) # NOT _under!
1296# _update_all(geodesic)
1297 _GeodesicLineExact.__init__(self, geodesic, lat1, lon1, azi1, caps, 0, **name)
1300def _Astroid(x, y):
1301 '''(INTERNAL) Solve M{k^4 + 2 * k^3 - (x^2 + y^2 - 1)
1302 * k^2 - (2 * k + 1) * y^2 = 0} for positive root k.
1303 '''
1304 p = x**2
1305 q = y**2
1306 r = fsumf_(_1_0, q, p, _N_2_0)
1307 if r > 0 or q:
1308 # avoid possible division by zero when r = 0
1309 # by multiplying s and t by r^3 and r, resp.
1310 S = p * q / _4_0 # S = r^3 * s
1311 if r:
1312 r = r / _6_0 # /= chokes PyChecker
1313 r3 = r**3
1314 T3 = r3 + S
1315 # discriminant of the quadratic equation for T3 is
1316 # zero on the evolute curve p^(1/3) + q^(1/3) = 1
1317 d = (r3 + T3) * S
1318 if d < 0:
1319 # T is complex, but u is defined for a real result
1320 a = atan2(sqrt(-d), -T3) / _3_0
1321 # There are 3 possible cube roots, choose the one which
1322 # avoids cancellation. Note d < 0 implies that r < 0.
1323 u = (cos(a) * _2_0 + _1_0) * r
1324 else:
1325 # pick the sign on the sqrt to maximize abs(T3) to
1326 # minimize loss of precision due to cancellation.
1327 if d:
1328 T3 += _copysign(sqrt(d), T3) # T3 = (r * t)^3
1329 # _cbrt always returns the real root, _cbrt(-8) = -2
1330 u = _cbrt(T3) # T = r * t
1331 if u: # T can be zero; but then r2 / T -> 0
1332 u += r**2 / u
1333 u += r
1334 elif S: # d == T3**2 == S**2: sqrt(d) == abs(S) == abs(T3)
1335 u = _cbrt(S * _2_0) # == T3 + _copysign(abs(S), T3)
1336 else:
1337 u = _0_0
1338 v = _hypot(u, y) # sqrt(u**2 + q)
1339 # avoid loss of accuracy when u < 0
1340 u = (q / (v - u)) if u < 0 else (v + u)
1341 w = (u - q) / (v + v) # positive?
1342 # rearrange expression for k to avoid loss of accuracy due to
1343 # subtraction, division by 0 impossible because u > 0, w >= 0
1344 k = u / (sqrt(w**2 + u) + w) # guaranteed positive
1346 else: # q == 0 && r <= 0
1347 # y = 0 with |x| <= 1. Handle this case directly, for
1348 # y small, positive root is k = abs(y) / sqrt(1 - x^2)
1349 k = _0_0
1351 return k
1354def _C4coeffs(nC4): # in .geodesicx.__main__
1355 '''(INTERNAL) Get the C{C4_24}, C{_27} or C{_30} series coefficients.
1356 '''
1357 try: # from pygeodesy.geodesicx._C4_xx import _coeffs_xx as _coeffs
1358 _C4_xx = _DOT_(_MODS.geodesicx.__name__, _UNDER_('_C4', nC4))
1359 _coeffs = _MODS.getattr(_C4_xx, _UNDER_('_coeffs', nC4))
1360 except (AttributeError, ImportError, TypeError) as x:
1361 raise GeodesicError(nC4=nC4, cause=x)
1362 n = _xnC4(nC4=nC4)
1363 if len(_coeffs) != n: # double check
1364 raise GeodesicError(_coeffs=len(_coeffs), _xnC4=n, nC4=nC4)
1365 return _coeffs
1368__all__ += _ALL_DOCS(GeodesicExact, GeodesicLineExact)
1370# **) MIT License
1371#
1372# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1373#
1374# Permission is hereby granted, free of charge, to any person obtaining a
1375# copy of this software and associated documentation files (the "Software"),
1376# to deal in the Software without restriction, including without limitation
1377# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1378# and/or sell copies of the Software, and to permit persons to whom the
1379# Software is furnished to do so, subject to the following conditions:
1380#
1381# The above copyright notice and this permission notice shall be included
1382# in all copies or substantial portions of the Software.
1383#
1384# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1385# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1386# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1387# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1388# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1389# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1390# OTHER DEALINGS IN THE SOFTWARE.