Coverage for pygeodesy/sphericalBase.py: 92%
211 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
2# -*- coding: utf-8 -*-
4u'''(INTERNAL) Private spherical base classes C{CartesianSphericalBase} and
5C{LatLonSphericalBase} for L{sphericalNvector} and L{sphericalTrigonometry}.
7A pure Python implementation of geodetic (lat-/longitude) functions,
8transcoded in part from JavaScript originals by I{(C) Chris Veness 2011-2016}
9and published under the same MIT Licence**, see
10U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.
11'''
12# make sure int/int division yields float quotient, see .basics
13from __future__ import division as _; del _ # PYCHOK semicolon
15from pygeodesy.basics import isbool, isinstanceof, map1
16from pygeodesy.cartesianBase import Bearing2Tuple, CartesianBase
17from pygeodesy.constants import EPS, PI, PI2, PI_2, R_M, R_MA, \
18 _umod_360, isnear0, isnon0, _0_0, \
19 _0_5, _1_0, _180_0
20from pygeodesy.datums import Datums, _spherical_datum
21from pygeodesy.errors import IntersectionError, _ValueError, _xError
22from pygeodesy.fmath import favg, fdot, hypot, sqrt_a
23from pygeodesy.interns import NN, _COMMA_, _concentric_, _datum_, \
24 _distant_, _exceed_PI_radians_, _name_, \
25 _near_, _radius_, _too_
26from pygeodesy.latlonBase import LatLonBase, _trilaterate5 # PYCHOK passed
27from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
28# from pygeodesy.namedTuples import Bearing2Tuple # from .cartesianBase
29from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs
30from pygeodesy.props import deprecated_method, property_doc_, \
31 property_RO, _update_all
32# from pygeodesy.streprs import Fmt, _xattrs # from .nvectorBase
33from pygeodesy.units import Bearing, Bearing_, Height, Radians_, \
34 Radius, Radius_, Scalar_
35from pygeodesy.utily import acos1, atan2b, atan2d, degrees90, degrees180, \
36 sincos2, sincos2d, tanPI_2_2, wrap360, wrapPI
38from math import cos, fabs, log, sin, sqrt
40__all__ = _ALL_LAZY.sphericalBase
41__version__ = '23.04.10'
44def _angular(distance, radius, low=EPS): # PYCHOK in .spherical*
45 '''(INTERNAL) Return the angular distance in C{radians}.
47 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}.
48 '''
49 r = _1_0 if radius is None else Radius_(radius=radius)
50 return Radians_(distance / r, low=low)
53def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry
54 '''(INTERNAL) Convert radii to radians.
55 '''
56 r1 = Radius_(rad1=rad1)
57 r2 = Radius_(rad2=rad2)
58 if radius is not None: # convert radii to radians
59 r1 = _angular(r1, radius)
60 r2 = _angular(r2, radius)
62 x = r1 < r2
63 if x:
64 r1, r2 = r2, r1
65 if r1 > PI:
66 raise IntersectionError(rad1=rad1, rad2=rad2,
67 txt=_exceed_PI_radians_)
68 return r1, r2, x
71class CartesianSphericalBase(CartesianBase):
72 '''(INTERNAL) Base class for spherical C{Cartesian}s.
73 '''
74 _datum = Datums.Sphere # L{Datum}
76 def intersections2(self, rad1, other, rad2, radius=R_M):
77 '''Compute the intersection points of two circles each defined
78 by a center point and a radius.
80 @arg rad1: Radius of the this circle (C{meter} or C{radians},
81 see B{C{radius}}).
82 @arg other: Center of the other circle (C{Cartesian}).
83 @arg rad2: Radius of the other circle (C{meter} or C{radians},
84 see B{C{radius}}).
85 @kwarg radius: Mean earth radius (C{meter} or C{None} if both
86 B{C{rad1}} and B{C{rad2}} are given in C{radians}).
88 @return: 2-Tuple of the intersection points, each C{Cartesian}.
89 For abutting circles, the intersection points are the
90 same C{Cartesian} instance, aka the I{radical center}.
92 @raise IntersectionError: Concentric, antipodal, invalid or
93 non-intersecting circles.
95 @raise TypeError: If B{C{other}} is not C{Cartesian}.
97 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}.
99 @see: U{Calculating intersection of two Circles
100 <https://GIS.StackExchange.com/questions/48937/
101 calculating-intersection-of-two-circles>} and method
102 or function C{trilaterate3d2}.
103 '''
104 x1, x2 = self, self.others(other)
105 r1, r2, x = _rads3(rad1, rad2, radius)
106 if x:
107 x1, x2 = x2, x1
108 try:
109 n, q = x1.cross(x2), x1.dot(x2)
110 n2, q1 = n.length2, (_1_0 - q**2)
111 if n2 < EPS or isnear0(q1):
112 raise ValueError(_near_(_concentric_))
113 c1, c2 = cos(r1), cos(r2)
114 x0 = x1.times((c1 - q * c2) / q1).plus(
115 x2.times((c2 - q * c1) / q1))
116 n1 = _1_0 - x0.length2
117 if n1 < EPS:
118 raise ValueError(_too_(_distant_))
119 except ValueError as x:
120 raise IntersectionError(center=self, rad1=rad1,
121 other=other, rad2=rad2, cause=x)
122 n = n.times(sqrt(n1 / n2))
123 if n.length > EPS:
124 x1 = x0.plus(n)
125 x2 = x0.minus(n)
126 else: # abutting circles
127 x1 = x2 = x0
129 return (_xattrs(x1, self, _datum_, _name_),
130 _xattrs(x2, self, _datum_, _name_))
133class LatLonSphericalBase(LatLonBase):
134 '''(INTERNAL) Base class for spherical C{LatLon}s.
135 '''
136 _datum = Datums.Sphere # spherical L{Datum}
138 def __init__(self, lat, lon, height=0, datum=None, name=NN):
139 '''Create a spherical C{LatLon} point frome the given
140 lat-, longitude and height on the given datum.
142 @arg lat: Latitude (C{degrees} or DMS C{[N|S]}).
143 @arg lon: Longitude (C{degrees} or DMS C{str[E|W]}).
144 @kwarg height: Optional elevation (C{meter}, the same units
145 as the datum's half-axes).
146 @kwarg datum: Optional, spherical datum to use (L{Datum},
147 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple})
148 or C{scalar} earth radius).
149 @kwarg name: Optional name (string).
151 @raise TypeError: If B{C{datum}} invalid or not
152 not spherical.
154 @example:
156 >>> p = LatLon(51.4778, -0.0016) # height=0, datum=Datums.WGS84
157 '''
158 LatLonBase.__init__(self, lat, lon, height=height, name=name)
159 if datum not in (None, self.datum):
160 self.datum = datum
162 def bearingTo2(self, other, wrap=False, raiser=False):
163 '''Return the initial and final bearing (forward and reverse
164 azimuth) from this to an other point.
166 @arg other: The other point (C{LatLon}).
167 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
168 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}).
170 @return: A L{Bearing2Tuple}C{(initial, final)}.
172 @raise TypeError: The B{C{other}} point is not spherical.
174 @see: Methods C{initialBearingTo} and C{finalBearingTo}.
175 '''
176 # .initialBearingTo is inside .-Nvector and .-Trigonometry
177 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo
178 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser)
179 return Bearing2Tuple(i, f, name=self.name)
181 @property_doc_(''' this point's datum (L{Datum}).''')
182 def datum(self):
183 '''Get this point's datum (L{Datum}).
184 '''
185 return self._datum
187 @datum.setter # PYCHOK setter!
188 def datum(self, datum):
189 '''Set this point's datum I{without conversion} (L{Datum}, L{Ellipsoid},
190 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} spherical earth radius).
192 @raise TypeError: If B{C{datum}} invalid or not not spherical.
193 '''
194 d = _spherical_datum(datum, name=self.name, raiser=_datum_)
195 if self._datum != d:
196 _update_all(self)
197 self._datum = d
199 def finalBearingTo(self, other, wrap=False, raiser=False):
200 '''Return the final bearing (reverse azimuth) from this to
201 an other point.
203 @arg other: The other point (spherical C{LatLon}).
204 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
205 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}).
207 @return: Final bearing (compass C{degrees360}).
209 @raise TypeError: The B{C{other}} point is not spherical.
211 @example:
213 >>> p = LatLon(52.205, 0.119)
214 >>> q = LatLon(48.857, 2.351)
215 >>> b = p.finalBearingTo(q) # 157.9
216 '''
217 self.others(other)
219 # final bearing is the reverse of the other, initial one;
220 # .initialBearingTo is inside .-Nvector and .-Trigonometry
221 b = other.initialBearingTo(self, wrap=wrap, raiser=raiser)
222 return _umod_360(b + _180_0)
224 def intersecant2(self, circle, point, bearing, radius=R_M, exact=False,
225 height=None, wrap=True):
226 '''Compute the intersections of a circle and a line.
228 @arg circle: Radius of the circle centered at this location
229 (C{meter}, same units as B{C{radius}}) or a point
230 on the circle (this C{LatLon}).
231 @arg point: An other point in- or outside the circle (this C{LatLon}).
232 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360})
233 or a second point on the line (this C{LatLon}).
234 @kwarg radius: Mean earth radius (C{meter}, conventionally).
235 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
236 destination and distance, if C{False} use the basic
237 rhumb methods (C{bool}) or if C{None} use the I{great
238 circle} methods.
239 @kwarg height: Optional height for the intersection points (C{meter},
240 conventionally) or C{None}.
241 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
243 @return: 2-Tuple of the intersection points (representing a chord),
244 each an instance of this class. For a tangent line, each
245 point C{is} this very instance.
247 @raise IntersectionError: The circle and line do not intersect.
249 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}
250 or B{C{bearing}} invalid.
252 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
253 B{C{exact}} or B{C{height}}.
254 '''
255 p = self.others(point=point)
256 try:
257 return _intersecant2(self, circle, p, bearing, radius=radius, exact=exact,
258 height=height, wrap=wrap)
259 except (TypeError, ValueError) as x:
260 raise _xError(x, center=self, circle=circle, point=point, bearing=bearing, exact=exact)
262 def maxLat(self, bearing):
263 '''Return the maximum latitude reached when travelling
264 on a great circle on given bearing from this point
265 based on Clairaut's formula.
267 The maximum latitude is independent of longitude
268 and the same for all points on a given latitude.
270 Negate the result for the minimum latitude (on the
271 Southern hemisphere).
273 @arg bearing: Initial bearing (compass C{degrees360}).
275 @return: Maximum latitude (C{degrees90}).
277 @raise ValueError: Invalid B{C{bearing}}.
278 '''
279 m = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi)))
280 return degrees90(m)
282 def minLat(self, bearing):
283 '''Return the minimum latitude reached when travelling
284 on a great circle on given bearing from this point.
286 @arg bearing: Initial bearing (compass C{degrees360}).
288 @return: Minimum latitude (C{degrees90}).
290 @see: Method L{maxLat} for more details.
292 @raise ValueError: Invalid B{C{bearing}}.
293 '''
294 return -self.maxLat(bearing)
296 def parse(self, strllh, height=0, sep=_COMMA_, name=NN):
297 '''Parse a string representing a similar, spherical C{LatLon}
298 point, consisting of C{"lat, lon[, height]"}.
300 @arg strllh: Lat, lon and optional height (C{str}),
301 see function L{pygeodesy.parse3llh}.
302 @kwarg height: Optional, default height (C{meter}).
303 @kwarg sep: Optional separator (C{str}).
304 @kwarg name: Optional instance name (C{str}),
305 overriding this name.
307 @return: The similar point (spherical C{LatLon}).
309 @raise ParseError: Invalid B{C{strllh}}.
310 '''
311 t = _MODS.dms.parse3llh(strllh, height=height, sep=sep)
312 r = self.classof(*t)
313 if name:
314 r.rename(name)
315 return r
317 @property_RO
318 def _radius(self):
319 '''(INTERNAL) Get this sphere's radius.
320 '''
321 return self.datum.ellipsoid.equatoradius
323 def _rhumb3(self, other, r=False):
324 '''(INTERNAL) Rhumb_ helper function.
326 @arg other: The other point (spherical C{LatLon}).
327 '''
328 self.others(other)
330 a1, b1 = self.philam
331 a2, b2 = other.philam
332 # if |db| > 180 take shorter rhumb
333 # line across the anti-meridian
334 db = wrapPI(b2 - b1)
335 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1))
336 da = a2 - a1
337 if r:
338 # on Mercator projection, longitude distances shrink
339 # by latitude; the 'stretch factor' q becomes ill-
340 # conditioned along E-W line (0/0); use an empirical
341 # tolerance to avoid it
342 q = (da / dp) if fabs(dp) > EPS else cos(self.phi)
343 da = hypot(da, q * db) # angular distance radians
344 return da, db, dp
346 def rhumbAzimuthTo(self, other, radius=R_M, exact=False):
347 '''Return the azimuth (bearing) of a rhumb line (loxodrome)
348 between this and an other (spherical) point.
350 @arg other: The other point (spherical C{LatLon}).
351 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
352 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).
353 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default
354 C{False} for backward compatibility.
356 @return: Rhumb line azimuth (compass C{degrees180}).
358 @raise TypeError: The B{C{other}} point is incompatible or
359 B{C{radius}} is invalid.
361 @example:
363 >>> p = LatLon(51.127, 1.338)
364 >>> q = LatLon(50.964, 1.853)
365 >>> b = p.rhumbBearingTo(q) # 116.7
366 '''
367 if exact: # use series, always
368 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False, radius=radius)
369 else:
370 _, db, dp = self._rhumb3(other)
371 z = atan2d(db, dp) # see .rhumbx.Rhumb.Inverse
372 return z
374 @deprecated_method
375 def rhumbBearingTo(self, other):
376 '''DEPRECATED, use method C{.rhumbAzimuthTo}.'''
377 return wrap360(self.rhumbAzimuthTo(other)) # [0..360)
379 def rhumbDestination(self, distance, bearing, radius=R_M, height=None, exact=False):
380 '''Return the destination point having travelled the given distance
381 from this point along a rhumb line (loxodrome) at the given bearing.
383 @arg distance: Distance travelled (C{meter}, same units as B{C{radius}}),
384 may be negative if C{B{exact}=True}.
385 @arg bearing: Bearing (azimuth) at this point (compass C{degrees360}).
386 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
387 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if
388 C{B{exact}=True}.
389 @kwarg height: Optional height, overriding the default height
390 (C{meter}, same unit as B{C{radius}}).
391 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default
392 C{False} for backward compatibility.
394 @return: The destination point (spherical C{LatLon}).
396 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
397 B{C{radius}} or B{C{height}}.
399 @example:
401 >>> p = LatLon(51.127, 1.338)
402 >>> q = p.rhumbDestination(40300, 116.7) # 50.9642°N, 001.8530°E
403 '''
404 if exact: # use series, always
405 r = LatLonBase.rhumbDestination(self, distance, bearing, exact=False,
406 radius=radius, height=height)
407 else: # radius=None from .rhumbMidpointTo
408 if radius in (None, self._radius):
409 d, r = self.datum, radius
410 else:
411 d = _spherical_datum(radius, raiser=_radius_) # spherical only
412 r = d.ellipsoid.equatoradius
413 r = _angular(distance, r, low=_0_0) # distance=0 from .rhumbMidpointTo
415 a1, b1 = self.philam
416 sb, cb = sincos2(Bearing_(bearing))
418 da = r * cb
419 a2 = a1 + da
420 # normalize latitude if past pole
421 if a2 > PI_2:
422 a2 = PI - a2
423 elif a2 < -PI_2:
424 a2 = -PI - a2
426 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1))
427 # q becomes ill-conditioned on E-W course 0/0
428 q = (da / dp) if fabs(dp) > EPS else cos(a1)
429 b2 = (b1 + r * sb / q) if fabs(q) > EPS else b1
431 h = self.height if height is None else Height(height)
432 r = self.classof(degrees90(a2), degrees180(b2), datum=d, height=h)
433 return r
435 def rhumbDistanceTo(self, other, radius=R_M, exact=False):
436 '''Return the distance from this to an other point along
437 a rhumb line (loxodrome).
439 @arg other: The other point (spherical C{LatLon}).
440 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
441 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if
442 C{B{exact}=True}.
443 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default
444 C{False} for backward compatibility.
446 @return: Distance (C{meter}, the same units as B{C{radius}}
447 or C{radians} if B{C{radius}} is C{None}).
449 @raise TypeError: The B{C{other}} point is incompatible.
451 @raise ValueError: Invalid B{C{radius}}.
453 @example:
455 >>> p = LatLon(51.127, 1.338)
456 >>> q = LatLon(50.964, 1.853)
457 >>> d = p.rhumbDistanceTo(q) # 403100
458 '''
459 if exact: # use series, always
460 r = LatLonBase.rhumbDistanceTo(self, other, exact=False, radius=radius)
461 if radius is None: # angular distance in radians
462 r = r / self._radius # /= chokes PyChecker
463 else:
464 # see <https://www.EdWilliams.org/avform.htm#Rhumb>
465 r, _, _ = self._rhumb3(other, r=True)
466 if radius is not None:
467 r *= Radius(radius)
468 return r
470 def rhumbMidpointTo(self, other, height=None, radius=R_M,
471 exact=False, fraction=_0_5):
472 '''Return the (loxodromic) midpoint on the rhumb line between
473 this and an other point.
475 @arg other: The other point (spherical LatLon).
476 @kwarg height: Optional height, overriding the mean height
477 (C{meter}).
478 @kwarg radius: Optional mean earth radius (C{meter}),
479 overriding the default C{R_M}.
480 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
481 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).
482 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default
483 C{False} for backward compatibility.
484 @kwarg fraction: Midpoint location from this point (C{scalar}),
485 may negative if C{B{exact}=True}.
487 @return: The (mid)point at the given B{C{fraction}} along
488 the rhumb line (spherical C{LatLon}).
490 @raise TypeError: The B{C{other}} point is incompatible.
492 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}
494 @example:
496 >>> p = LatLon(51.127, 1.338)
497 >>> q = LatLon(50.964, 1.853)
498 >>> m = p.rhumb_midpointTo(q)
499 >>> m.toStr() # '51.0455°N, 001.5957°E'
500 '''
501 if exact: # use series, always
502 r = LatLonBase.rhumbMidpointTo(self, other, exact=False,
503 radius=radius, height=height, fraction=fraction)
504 elif fraction is not _0_5:
505 f = Scalar_(fraction=fraction) # low=_0_0
506 r, db, dp = self._rhumb3(other, r=True) # radians
507 z = atan2b(db, dp)
508 h = self._havg(other, f=f) if height is None else height
509 r = self.rhumbDestination(r * f, z, radius=None, height=h)
511 else: # for backward compatibility
512 self.others(other)
513 # see <https://MathForum.org/library/drmath/view/51822.html>
514 a1, b1 = self.philam
515 a2, b2 = other.philam
516 if fabs(b2 - b1) > PI:
517 b1 += PI2 # crossing anti-meridian
519 a3 = favg(a1, a2)
520 b3 = favg(b1, b2)
522 f1 = tanPI_2_2(a1)
523 if isnon0(f1):
524 f2 = tanPI_2_2(a2)
525 f = f2 / f1
526 if isnon0(f):
527 f = log(f)
528 if isnon0(f):
529 f3 = tanPI_2_2(a3)
530 b3 = fdot(map1(log, f1, f2, f3),
531 -b2, b1, b2 - b1) / f
533 d = self.datum if radius in (None, self._radius) else \
534 _spherical_datum(radius, name=self.name, raiser=_radius_)
535 h = self._havg(other) if height is None else Height(height)
536 r = self.classof(degrees90(a3), degrees180(b3), datum=d, height=h)
537 return r
539 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature
540 '''Convert this point to C{Nvector} components, I{including
541 height}.
543 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}
544 keyword arguments, ignored if
545 C{B{Nvector} is None}.
547 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}
548 if B{C{Nvector}} is C{None}.
550 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.
551 '''
552 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)
554 def toWm(self, radius=R_MA):
555 '''Convert this point to a I{WM} coordinate.
557 @kwarg radius: Optional earth radius (C{meter}).
559 @return: The WM coordinate (L{Wm}).
561 @see: Function L{pygeodesy.toWm} in module L{webmercator} for details.
562 '''
563 return _MODS.webmercator.toWm(self, radius=radius)
566def _intersecant2(c, r, p, b, radius=R_M, exact=False,
567 height=None, wrap=False):
568 # (INTERNAL) Intersect a circle and bearing, see L{intersecant2}
569 # above, separated to allow callers to embellish any exceptions
571 if not isinstanceof(r, c.__class__, p.__class__):
572 r = Radius_(circle=r)
573 elif exact is None:
574 r = c.distanceTo(r, radius=radius, wrap=wrap)
575 elif isbool(exact):
576 r = c.rhumbDistanceTo(r, radius=radius, exact=exact)
577 else:
578 raise _ValueError(exact=exact)
580 if not isinstanceof(b, c.__class__, p.__class__):
581 b = Bearing(b)
582 elif exact is None:
583 b = p.initialBearingTo(b, wrap=wrap)
584 else:
585 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact)
587 if exact is None:
588 d = p.distanceTo(c, radius=radius, wrap=wrap)
589 a = p.initialBearingTo(c, wrap=wrap) if d > EPS else b
590 else:
591 d = p.rhumbDistanceTo(c, radius=radius, exact=exact)
592 a = wrap360(p.rhumbAzimuthTo( c, radius=radius, exact=exact)) if d > EPS else b
594 if d > EPS:
595 s, c = sincos2d(b - a)
596 s = sqrt_a(r, fabs(s * d))
597 if s > r:
598 raise IntersectionError(_too_(Fmt.distant(s)))
599 elif (r - s) < EPS:
600 return p, p # tangent
601 c *= d
602 else: # coincindent
603 s, c = r, 0
604 a = b + _180_0
605 if exact is None:
606 b = p.destination(s + c, b, radius=radius, height=height)
607 a = p.destination(s - c, a, radius=radius, height=height)
608 else:
609 b = p.rhumbDestination(s + c, b, radius=radius, height=height, exact=exact)
610 a = p.rhumbDestination(s - c, a, radius=radius, height=height, exact=exact)
611 return b, a # in bearing direction first
614__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase)
616# **) MIT License
617#
618# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
619#
620# Permission is hereby granted, free of charge, to any person obtaining a
621# copy of this software and associated documentation files (the "Software"),
622# to deal in the Software without restriction, including without limitation
623# the rights to use, copy, modify, merge, publish, distribute, sublicense,
624# and/or sell copies of the Software, and to permit persons to whom the
625# Software is furnished to do so, subject to the following conditions:
626#
627# The above copyright notice and this permission notice shall be included
628# in all copies or substantial portions of the Software.
629#
630# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
631# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
632# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
633# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
634# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
635# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
636# OTHER DEALINGS IN THE SOFTWARE.