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