Coverage for pygeodesy/geodesicw.py: 89%
222 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
2# -*- coding: utf-8 -*-
4u'''Wrapper around Python classes C{geodesic.Geodesic} and C{geodesicline.GeodesicLine} from
5I{Karney}'s Python package U{geographiclib<https://PyPI.org/project/geographiclib>}, provided
6that package is installed.
8The I{wrapped} class methods return a L{GDict} instance offering access to the C{dict} items
9either by C{key} or by C{attribute} name.
11With env variable C{PYGEODESY_GEOGRAPHICLIB} left undefined or set to C{"2"}, this module,
12L{pygeodesy.geodesicx} and L{pygeodesy.karney} will use U{GeographicLib 2.0
13<https://GeographicLib.SourceForge.io/C++/doc/>} transcoding, otherwise C{1.52} or older.
14'''
16from pygeodesy.basics import _copysign, _xinstanceof
17from pygeodesy.constants import EPS, NAN, _EPSqrt as _TOL, _0_5
18from pygeodesy.datums import _earth_datum, _WGS84, _EWGS84
19# from pygeodesy.dms import F_D # from .latlonBase
20# from pygeodesy.ellipsoids import _EWGS84 # from .datums
21from pygeodesy.errors import IntersectionError, GeodesicError
22from pygeodesy.fsums import Fsum, Fmt, unstr
23from pygeodesy.internals import _dunder_nameof, _under
24from pygeodesy.interns import NN, _DOT_, _SPACE_, _to_, _too_
25from pygeodesy.karney import _atan2d, Caps, Direct9Tuple, GDict, \
26 _kWrapped, Inverse10Tuple
27from pygeodesy.latlonBase import LatLonBase as _LLB, F_D, Radius_
28from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
29from pygeodesy.named import callername, classname, _name1__, _name2__
30from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple
31from pygeodesy.props import Property, Property_RO, property_RO
32# from pygeodesy.streprs import Fmt, unstr # from .fsums
33# from pygeodesy.units import Radius_ # from .latlonBase
34from pygeodesy.utily import _unrollon, _Wrap, wrap360, fabs # PYCHOK used!
36from contextlib import contextmanager
37# from math import fabs # from .utily
39__all__ = _ALL_LAZY.geodesicw
40__version__ = '24.05.24'
42_plumb_ = 'plumb'
43_TRIPS = 65
46class _gWrapped(_kWrapped):
47 ''''(INTERNAL) Wrapper for some of I{Karney}'s U{geographiclib
48 <https://PyPI.org/project/geographiclib>} classes.
49 '''
51 @Property_RO # MCCABE 24
52 def Geodesic(self):
53 '''Get the I{wrapped} C{geodesic.Geodesic} class from I{Karney}'s Python
54 U{geographiclib<https://GitHub.com/geographiclib/geographiclib-python>},
55 provided the latter is installed.
56 '''
57 _Geodesic = self.geographiclib.Geodesic
58 # assert Caps._STD == _Geodesic.STANDARD
60 class Geodesic(_Geodesic):
61 '''I{Wrapper} for I{Karney}'s Python U{geodesic.Geodesic
62 <https://PyPI.org/project/geographiclib>} class.
63 '''
64 _datum = _WGS84
65 _debug = 0 # like .geodesicx.bases._GeodesicBase
66 LINE_OFF = 0 # in .azimuthal._GnomonicBase and .css.CassiniSoldner
67 _name = NN
69 def __init__(self, a_ellipsoid=_EWGS84, f=None, **name): # PYCHOK signature
70 '''New I{wrapped} C{geodesic.Geodesic} instance.
72 @arg a_ellipsoid: The equatorial radius I{a} (C{meter}, conventionally),
73 an ellipsoid (L{Ellipsoid}) or a datum (L{Datum}).
74 @arg f: The ellipsoid's flattening (C{scalar}), ignored if B{C{a_ellipsoid})
75 is not C{meter}.
76 @kwarg name: Optional C{B{name}=NN} (C{str}).
77 '''
78 _earth_datum(self, a_ellipsoid, f=f, **name) # raiser=NN
79 E = self.ellipsoid
80 with _wargs(self, *E.a_f, **name) as args:
81 _Geodesic.__init__(self, *args)
82 if name:
83 self._name, _ = _name2__(name, _or_nameof=E)
85 def ArcDirect(self, lat1, lon1, azi1, a12, outmask=Caps._STD):
86 '''Return the C{_Geodesic.ArcDirect} result as L{GDict}.
87 '''
88 with _wargs(self, lat1, lon1, azi1, a12, outmask) as args:
89 d = _Geodesic.ArcDirect(self, *args)
90 return GDict(d)
92 def ArcDirectLine(self, lat1, lon1, azi1, a12, caps=Caps._STD_LINE, **name):
93 '''Return the C{_Geodesic.ArcDirectLine} as I{wrapped} C{GeodesicLine}.
94 '''
95 return self._GenDirectLine(lat1, lon1, azi1, True, a12, caps, **name)
97 Area = _Geodesic.Polygon # like GeodesicExact.Area
99 @property_RO
100 def datum(self):
101 '''Get this geodesic's datum (C{Datum}).
102 '''
103 return self._datum
105 @Property
106 def debug(self):
107 '''Get the C{debug} option (C{bool}).
108 '''
109 return bool(self._debug)
111 @debug.setter # PYCHOK setter!
112 def debug(self, debug):
113 '''Set the C{debug} option (C{bool}) to include more
114 details in L{GDict} results.
115 '''
116 self._debug = Caps._DEBUG_ALL if debug else 0
118 def Direct(self, lat1, lon1, azi1, s12=0, outmask=Caps._STD):
119 '''Return the C{_Geodesic.Direct} result as L{GDict}.
120 '''
121 with _wargs(self, lat1, lon1, azi1, s12, outmask) as args:
122 d = _Geodesic.Direct(self, *args)
123 return GDict(d)
125 def Direct3(self, lat1, lon1, azi1, s12): # PYCHOK outmask
126 '''Return the destination lat, lon and reverse azimuth
127 in C{degrees} as L{Destination3Tuple}.
128 '''
129 d = self.Direct(lat1, lon1, azi1, s12, outmask=Caps._DIRECT3)
130 return Destination3Tuple(d.lat2, d.lon2, d.azi2)
132 def _DirectLine(self, ll1, azi12, s12=0, **caps_name):
133 '''(INTERNAL) Short-cut version.
134 '''
135 return self.DirectLine(ll1.lat, ll1.lon, azi12, s12, **caps_name)
137 def DirectLine(self, lat1, lon1, azi1, s12, caps=Caps._STD_LINE, **name):
138 '''Return the C{_Geodesic.DirectLine} as I{wrapped} C{GeodesicLine}.
139 '''
140 return self._GenDirectLine(lat1, lon1, azi1, False, s12, caps, **name)
142 @Property_RO
143 def ellipsoid(self):
144 '''Get this geodesic's ellipsoid (C{Ellipsoid}).
145 '''
146 return self.datum.ellipsoid
148 @property_RO
149 def f1(self): # in .css.CassiniSoldner.reset
150 '''Get the geodesic's ellipsoid's I{1 - flattening} (C{float}).
151 '''
152 return getattr(self, _under(Geodesic.f1.name), self.ellipsoid.f1)
154 def _GDictDirect(self, lat, lon, azi, arcmode, s12_a12, outmask=Caps._STD):
155 '''(INTERNAL) Get C{_Geodesic._GenDirect} result as C{GDict}.
156 '''
157 with _wargs(self, lat, lon, azi, arcmode, s12_a12, outmask) as args:
158 t = _Geodesic._GenDirect(self, *args)
159 return Direct9Tuple(t).toGDict() # *t
161 def _GDictInverse(self, lat1, lon1, lat2, lon2, outmask=Caps._STD):
162 '''(INTERNAL) Get C{_Geodesic._GenInverse} result as L{Inverse10Tuple}.
163 '''
164 with _wargs(self, lat1, lon1, lat2, lon2, outmask) as args:
165 t = _Geodesic._GenInverse(self, *args)
166 return Inverse10Tuple(t).toGDict(lon1=lon1, lon2=lon2) # *t
168 def _GenDirectLine(self, lat1, lon1, azi1, arcmode, s12_a12, *caps, **name):
169 '''(INTERNAL) Invoked by C{_Geodesic.DirectLine} and C{-.ArcDirectLine},
170 returning the result as a I{wrapped} C{GeodesicLine}.
171 '''
172 with _wargs(self, lat1, lon1, azi1, arcmode, s12_a12, *caps, **name) as args:
173 t = _Geodesic._GenDirectLine(self, *args)
174 return self._Line13(t, **name)
176 def _Inverse(self, ll1, ll2, wrap, **outmask):
177 '''(INTERNAL) Short-cut version, see .ellipsoidalBaseDI.intersecant2.
178 '''
179 if wrap:
180 ll2 = _unrollon(ll1, _Wrap.point(ll2))
181 return self.Inverse(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **outmask)
183 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps._STD):
184 '''Return the C{_Geodesic.Inverse} result as L{GDict}.
185 '''
186 with _wargs(self, lat1, lon1, lat2, lon2, outmask) as args:
187 d = _Geodesic.Inverse(self, *args)
188 return GDict(d)
190 def Inverse1(self, lat1, lon1, lat2, lon2, wrap=False):
191 '''Return the non-negative, I{angular} distance in C{degrees}.
193 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
194 B{C{lat2}} and BC{lon2}} (C{bool}).
195 '''
196 # see .FrechetKarney.distance, .HausdorffKarney._distance
197 # and .HeightIDWkarney._distances
198 if wrap:
199 _, lat2, lon2 = _Wrap.latlon3(lat1, lat2, lon2, True) # _Geodesic.LONG_UNROLL
200 r = self.Inverse(lat1, lon1, lat2, lon2)
201 # XXX _Geodesic.DISTANCE needed for 'a12'?
202 return fabs(r.a12)
204 def Inverse3(self, lat1, lon1, lat2, lon2): # PYCHOK outmask
205 '''Return the distance in C{meter} and the forward and reverse
206 azimuths in C{degrees} as L{Distance3Tuple}.
207 '''
208 r = self.Inverse(lat1, lon1, lat2, lon2, outmask=Caps._INVERSE3)
209 return Distance3Tuple(r.s12, wrap360(r.azi1), wrap360(r.azi2))
211 def _InverseLine(self, ll1, ll2, wrap, **caps_name):
212 '''(INTERNAL) Short-cut version.
213 '''
214 if wrap:
215 ll2 = _unrollon(ll1, _Wrap.point(ll2))
216 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **caps_name)
218 def InverseLine(self, lat1, lon1, lat2, lon2, caps=Caps._STD_LINE, **name):
219 '''Return the C{_Geodesic.InverseLine} as I{wrapped} C{GeodesicLine}.
220 '''
221 with _wargs(self, lat1, lon1, lat2, lon2, caps, **name) as args:
222 t = _Geodesic.InverseLine(self, *args)
223 return self._Line13(t, **name)
225 def Line(self, lat1, lon1, azi1, caps=Caps._STD_LINE, **name):
226 '''Set up a I{wrapped} C{GeodesicLine} to compute several points
227 along a single, I{wrapped} (this) geodesic.
228 '''
229 return _wrapped.GeodesicLine(self, lat1, lon1, azi1, caps=caps, **name)
231 def _Line13(self, t, **name):
232 '''(INTERNAL) Wrap C{_GeodesicLine}, add distance and arc length
233 to reference point 3.
234 '''
235 gl = _wrapped.GeodesicLine(self, t.lat1, t.lon1, t.azi1, caps=t.caps,
236 salp1=t.salp1, calp1=t.calp1, **name)
237 gl.a13, gl.s13 = t.a13, t.s13
238 return gl
240 @property_RO
241 def name(self):
242 '''Get the name (C{str}).
243 '''
244 return self._name
246# Polygon = _Geodesic.Polygon
248 # Geodesic.ArcDirect.__doc__ = _Geodesic.ArcDirect.__doc__
249 # Geodesic.Direct.__doc__ = _Geodesic.Direct.__doc__
250 # Geodesic.Inverse.__doc__ = _Geodesic.Inverse.__doc__
251 # Geodesic.InverseLine.__doc__ = _Geodesic.InverseLinr.__doc__
252 # Geodesic.Line.__doc__ = _Geodesic.Line.__doc__
253 return Geodesic
255 @Property_RO # MCCABE 16
256 def GeodesicLine(self):
257 '''Get the I{wrapped} C{geodesicline.GeodesicLine} class from I{Karney}'s
258 Python U{geographiclib<https://GitHub.com/geographiclib/geographiclib-python>},
259 provided the latter is installed.
260 '''
261 _GeodesicLine = self.geographiclib.GeodesicLine
263 class GeodesicLine(_GeodesicLine):
264 '''I{Wrapper} for I{Karney}'s Python U{geodesicline.GeodesicLine
265 <https://PyPI.org/project/geographiclib>} class.
266 '''
267 _geodesic = None
268 _name = NN
270 def __init__(self, geodesic, lat1, lon1, azi1, **caps_name_): # salp1=NAN, calp1=NAN
271 '''New I{wrapped} C{geodesicline.GeodesicLine} instance.
273 @arg geodesic: A I{wrapped} C{Geodesic} instance.
274 @arg lat1: Latitude of the first points (C{degrees}).
275 @arg lon1: Longitude of the first points (C{degrees}).
276 @arg azi1: Azimuth at the first points (compass C{degrees360}).
277 @kwarg caps_name_: Optional keyword arguments C{B{caps}=Caps.STANDARD},
278 a bit-or'ed combination of L{Caps} values specifying the
279 capabilities the C{GeodesicLine} instance should possess,
280 an optional C{B{name}=NN} plus C{salp1=NAN} and C{calp1=NAN}
281 for I{INTERNAL} use.
282 '''
283 _xinstanceof(_wrapped.Geodesic, geodesic=geodesic)
284 with _wargs(self, geodesic, lat1, lon1, azi1, **caps_name_) as args:
285 name, caps_ = _name2__(caps_name_, _or_nameof=geodesic)
286 _GeodesicLine.__init__(self, *args, **caps_) # XXX avoid updates?
287 if name:
288 self._name = name
289 self._geodesic = geodesic
291 @Property_RO
292 def a1(self):
293 '''Get the I{equatorial arc} (C{degrees}), the arc length between
294 the northward equatorial crossing and point C{(lat1, lon1)}.
296 @see: U{EquatorialArc<https://GeographicLib.SourceForge.io/
297 C++/doc/classGeographicLib_1_1GeodesicLine.html>}
298 '''
299 try:
300 return _atan2d(self._ssig1, self._csig1)
301 except AttributeError:
302 return NAN # see .geodesicx.gxline._GeodesicLineExact
304 equatorarc = a1
306 def Arc(self):
307 '''Return the angular distance to point 3 (C{degrees} or C{NAN}).
308 '''
309 return self.a13
311 def ArcPosition(self, a12, outmask=Caps._STD):
312 '''Return the position at C{B{a12} degrees} on this line.
314 @arg a12: Angular distance from this line's first point
315 (C{degrees}).
317 @see: Method L{Position} for further details.
318 '''
319 with _wargs(self, a12, outmask) as args:
320 d = _GeodesicLine.ArcPosition(self, *args)
321 return GDict(d)
323 @Property_RO
324 def azi0(self): # see .css.CassiniSoldner.forward4
325 '''Get the I{equatorial azimuth} (C{degrees}), the azimuth of the
326 geodesic line as it crosses the equator in a northward direction.
328 @see: U{EquatorialAzimuth<https://GeographicLib.SourceForge.io/
329 C++/doc/classGeographicLib_1_1GeodesicLine.html>}
330 '''
331 try:
332 return _atan2d(self._salp0, self._calp0)
333 except AttributeError:
334 return NAN # see .geodesicx.gxline._GeodesicLineExact
336 equatorazimuth = azi0
338 def Distance(self):
339 '''Return the distance to reference point 3 (C{meter} or C{NAN}).
340 '''
341 return self.s13
343 @property_RO
344 def geodesic(self):
345 '''Get the I{wrapped} geodesic (L{Geodesic}).
346 '''
347 return self._geodesic
349 def Intersecant2(self, lat0, lon0, radius, tol=_TOL):
350 '''Compute the intersection(s) of this geodesic line and a circle.
352 @arg lat0: Latitude of the circle center (C{degrees}).
353 @arg lon0: Longitude of the circle center (C{degrees}).
354 @arg radius: Radius of the circle (C{meter}, conventionally).
355 @kwarg tol: Convergence tolerance (C{scalar}).
357 @return: 2-Tuple C{(P, Q)} with both intersections points (representing
358 a geodesic chord), each a L{GDict} from method L{Position} and
359 extended to 14 items C{lat1, lon1, azi1, lat2, lon2, azi2, a12,
360 s12, lat0, lon0, azi0, a02, s02, at} with the circle center
361 C{lat0}, C{lon0}, azimuth C{azi0} at the intersection, distance
362 C{a02} in C{degrees} and C{s02} in C{meter} along the geodesic
363 from the circle center to the intersection C{lat2, lon2} and
364 the angle C{at} between the geodesic and this line at the
365 intersection. The I{geodesic} azimuth at the intersection is
366 C{(at + azi2)}. If this line is tangential to the circle, both
367 intersections are the same L{GDict} instance.
369 @raise IntersectionError: The circle and this geodesic line do not
370 intersect.
372 @raise UnitError: Invalid B{C{radius}}.
373 '''
374 return _Intersecant2(self, lat0, lon0, radius, tol=tol)
376 def PlumbTo(self, lat0, lon0, est=None, tol=_TOL):
377 '''Compute the I{perpendicular} intersection of this geodesic line
378 with a geodesic from the given point.
380 @arg lat0: Latitude of the point (C{degrees}).
381 @arg lon0: Longitude of the point (C{degrees}).
382 @kwarg est: Optional, initial estimate for the distance C{s12} of
383 the intersection I{along} this geodesic line (C{meter}).
384 @kwarg tol: Convergence tolerance (C(meter)).
386 @return: The intersection point on this geodesic line, a L{GDict}
387 from method L{Position} extended to 14 items C{lat1, lon1,
388 azi1, lat2, lon2, azi2, a12, s12, lat0, lon0, azi0, a02,
389 s02, at} with C{a02} and C{s02} the distance in C{degrees}
390 and C{meter} from the given point C{lat0, lon0} to the
391 intersection C{lat2, lon2}, azimuth C{azi0} at the given
392 point and the (perpendicular) angle C{at} between the
393 geodesic and this line at the intersection point. The
394 geodesic azimuth at the intersection is C{(at + azi2)}.
395 See method L{Position} for further details.
397 @see: Methods C{Intersecant2}, C{Intersection} and C{Position}.
398 '''
399 return _PlumbTo(self, lat0, lon0, est=est, tol=tol)
401 def Position(self, s12, outmask=Caps._STD):
402 '''Return the position at distance C{B{s12} meter} on this line.
404 @arg s12: Distance from this line's first point (C{meter}).
405 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying
406 the quantities to be returned.
408 @return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2,
409 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1},
410 C{lon1}, C{azi1} and arc length C{a12} always included,
411 except when C{a12=NAN}.
412 '''
413 with _wargs(self, s12, outmask) as args:
414 d = _GeodesicLine.Position(self, *args)
415 return GDict(d)
417 # GeodesicLine.ArcPosition.__doc__ = _GeodesicLine.ArcPosition.__doc__
418 # GeodesicLine.Position.__doc__ = _GeodesicLine.Position.__doc__
419 return GeodesicLine
421 @Property_RO
422 def Geodesic_WGS84(self):
423 '''Get the I{wrapped} C{Geodesic(WGS84)} singleton, provided the
424 U{geographiclib<https://PyPI.org/project/geographiclib>} package
425 is installed, otherwise an C{ImportError}.
426 '''
427 return _EWGS84.geodesic
429_wrapped = _gWrapped() # PYCHOK singleton, .ellipsoids, .test/base.py
432def Geodesic(a_ellipsoid, f=None, **name):
433 '''Return a I{wrapped} C{geodesic.Geodesic} instance from I{Karney}'s
434 Python U{geographiclib<https://PyPI.org/project/geographiclib>},
435 provide the latter is installed, otherwise an C{ImportError}.
437 @arg a_ellipsoid: An ellipsoid (L{Ellipsoid}) or datum (L{Datum})
438 or the equatorial radius I{a} of the ellipsoid (C{meter}).
439 @arg f: The flattening of the ellipsoid (C{scalar}), ignored if
440 B{C{a_ellipsoid}}) is not specified as C{meter}.
441 @kwarg name: Optional ellipsoid C{B{name}=NN} (C{str}), ignored
442 like B{C{f}}.
443 '''
444 return _wrapped.Geodesic(a_ellipsoid, f=f, **name)
447def GeodesicLine(geodesic, lat1, lon1, azi1, caps=Caps._STD_LINE):
448 '''Return a I{wrapped} C{geodesicline.GeodesicLine} instance from I{Karney}'s
449 Python U{geographiclib<https://PyPI.org/project/geographiclib>}, provided
450 the latter is installed, otherwise an C{ImportError}.
452 @arg geodesic: A I{wrapped} L{Geodesic} instance.
453 @arg lat1: Latitude of the first points (C{degrees}).
454 @arg lon1: Longitude of the first points (C{degrees}).
455 @arg azi1: Azimuth at the first points (compass C{degrees360}).
456 @kwarg caps: Optional, bit-or'ed combination of L{Caps} values
457 specifying the capabilities the C{GeodesicLine}
458 instance should possess, i.e., which quantities can
459 be returned by calls to C{GeodesicLine.Position}
460 and C{GeodesicLine.ArcPosition}.
461 '''
462 return _wrapped.GeodesicLine(geodesic, lat1, lon1, azi1, caps=caps)
465def Geodesic_WGS84():
466 '''Get the I{wrapped} L{Geodesic}C{(WGS84)} singleton, provided
467 U{geographiclib<https://PyPI.org/project/geographiclib>} is
468 installed, otherwise an C{ImportError}.
469 '''
470 return _wrapped.Geodesic_WGS84
473class _wargs(object): # see also .formy._idllmn6, .latlonBase._toCartesian3, .vector2d._numpy
474 '''(INTERNAL) C{geographiclib} caller and exception mapper.
475 '''
476 @contextmanager # <https://www.Python.org/dev/peps/pep-0343/> Examples
477 def __call__(self, inst, *args, **kwds):
478 '''(INTERNAL) Yield C{tuple(B{args})} with any errors raised
479 as L{GeodesicError} embellished with all B{C{kwds}}.
480 '''
481 try:
482 yield args
483 except Exception as x:
484 u = _DOT_(classname(inst), callername(up=2, underOK=True))
485 raise GeodesicError(unstr(u, *args, **_name1__(kwds)), cause=x)
487_wargs = _wargs() # PYCHOK singleton
490def _Intersecant2(gl, lat0, lon0, radius, tol=_TOL, form=F_D): # MCCABE in LatLonEllipsoidalBaseDI.intersecant2, .geodesicx.gxline.Intersecant2
491 # (INTERNAL) Return the intersections of a circle at C{lat0, lon0}
492 # and a geodesic line as a 2-Tuple C{(P, Q)}, each a C{GDict}.
493 r = Radius_(radius)
494 n = _dunder_nameof(_Intersecant2)[1:]
495 _P = gl.Position
496 _I = gl.geodesic.Inverse
497 _a = fabs
499 def _R3(s):
500 # radius, intersection, etc. at distance C{s}
501 P = _P(s)
502 d = _I(lat0, lon0, P.lat2, P.lon2)
503 return _a(d.s12), P, d
505 def _bisect2(s, c, Rc, r, tol, _R3):
506 _s = Fsum(c).fsumf_
507 for i in range(_TRIPS):
508 b = _s(s)
509 Rb, P, d = _R3(b)
510 if Rb > r:
511 break
512 else: # b >>> s and c >>> s
513 raise ValueError(Fmt.no_convergence(b, s))
514 __2 = _0_5 # Rb > r > Rc
515 for i in range(_TRIPS): # 47-48
516 s = (b + c) * __2
517 R, P, d = _R3(s)
518 if Rb > R > r:
519 b, Rb = s, R
520 elif Rc < R < r:
521 c, Rc = s, R
522 t = _a(b - c)
523 if t < tol: # or _a(R - r) < tol:
524 break
525 else: # t = min(t, _a(R - r))
526 raise ValueError(Fmt.no_convergence(t, tol))
527 i += C.iteration # combine iterations
528 P.set_(lat0=lat0, lon0=lon0, azi0=d.azi1, iteration=i,
529 a02=d.a12, s02=d.s12, at=d.azi2 - P.azi2, name=n)
530 return P, s
532 # get the perpendicular intersection of 2 geodesics,
533 # one the plumb, pseudo-rhumb line to the other
534 C = _PlumbTo(gl, lat0, lon0, tol=tol)
535 try:
536 a = _a(C.s02) # distance between centers
537 if a < r:
538 c = C.s12 # distance along pseudo-rhumb line
539 h = _copysign(r, c) # past half chord length
540 P, p = _bisect2( h, c, a, r, tol, _R3)
541 Q, q = _bisect2(-h, c, a, r, tol, _R3)
542 if _a(p - q) < max(EPS, tol):
543 Q = P
544 elif a > r:
545 raise ValueError(_too_(Fmt.distant(a)))
546 else: # tangential
547 P = Q = C
548 except Exception as x:
549 t = _LLB(C.lat2, C.lon2).toStr(form=form)
550 t = _SPACE_(x, _plumb_, _to_, Fmt.PAREN(t))
551 raise IntersectionError(t, txt=None, cause=x)
553 return P, Q
556def _PlumbTo(gl, lat0, lon0, est=None, tol=_TOL):
557 # (INTERNAL) Return the I{perpendicular} intersection of
558 # a geodesic from C{(lat0, lon0)} and a geodesic (line).
559 pl = _MODS.rhumb.bases._PseudoRhumbLine(gl)
560 return pl.PlumbTo(lat0, lon0, exact=gl.geodesic,
561 est=est, tol=tol)
563# **) MIT License
564#
565# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
566#
567# Permission is hereby granted, free of charge, to any person obtaining a
568# copy of this software and associated documentation files (the "Software"),
569# to deal in the Software without restriction, including without limitation
570# the rights to use, copy, modify, merge, publish, distribute, sublicense,
571# and/or sell copies of the Software, and to permit persons to whom the
572# Software is furnished to do so, subject to the following conditions:
573#
574# The above copyright notice and this permission notice shall be included
575# in all copies or substantial portions of the Software.
576#
577# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
578# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
579# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
580# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
581# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
582# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
583# OTHER DEALINGS IN THE SOFTWARE.