Coverage for pygeodesy/sphericalBase.py: 91%
214 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-06-07 08:37 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-06-07 08:37 -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 CartesianBase, Bearing2Tuple
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_, Radians_, Radius, \
34 Radius_, Scalar_
35from pygeodesy.utily import acos1, atan2b, atan2d, degrees90, \
36 degrees180, sincos2, sincos2d, tanPI_2_2, \
37 _unrollon, wrap360, wrapPI
39from math import cos, fabs, log, sin, sqrt
41__all__ = _ALL_LAZY.sphericalBase
42__version__ = '23.06.05'
45def _angular(distance, radius, low=EPS): # PYCHOK in .spherical*
46 '''(INTERNAL) Return the angular distance in C{radians}.
48 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}.
49 '''
50 r = _1_0 if radius is None else Radius_(radius=radius)
51 return Radians_(distance / r, low=low)
54def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry
55 '''(INTERNAL) Convert radii to radians.
56 '''
57 r1 = Radius_(rad1=rad1)
58 r2 = Radius_(rad2=rad2)
59 if radius is not None: # convert radii to radians
60 r1 = _angular(r1, radius)
61 r2 = _angular(r2, radius)
63 x = r1 < r2
64 if x:
65 r1, r2 = r2, r1
66 if r1 > PI:
67 raise IntersectionError(rad1=rad1, rad2=rad2,
68 txt=_exceed_PI_radians_)
69 return r1, r2, x
72class CartesianSphericalBase(CartesianBase):
73 '''(INTERNAL) Base class for spherical C{Cartesian}s.
74 '''
75 _datum = Datums.Sphere # L{Datum}
77 def intersections2(self, rad1, other, rad2, radius=R_M):
78 '''Compute the intersection points of two circles each defined
79 by a center point and a radius.
81 @arg rad1: Radius of the this circle (C{meter} or C{radians},
82 see B{C{radius}}).
83 @arg other: Center of the other circle (C{Cartesian}).
84 @arg rad2: Radius of the other circle (C{meter} or C{radians},
85 see B{C{radius}}).
86 @kwarg radius: Mean earth radius (C{meter} or C{None} if both
87 B{C{rad1}} and B{C{rad2}} are given in C{radians}).
89 @return: 2-Tuple of the intersection points, each C{Cartesian}.
90 For abutting circles, the intersection points are the
91 same C{Cartesian} instance, aka the I{radical center}.
93 @raise IntersectionError: Concentric, antipodal, invalid or
94 non-intersecting circles.
96 @raise TypeError: If B{C{other}} is not C{Cartesian}.
98 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}.
100 @see: U{Calculating intersection of two Circles
101 <https://GIS.StackExchange.com/questions/48937/
102 calculating-intersection-of-two-circles>} and method
103 or function C{trilaterate3d2}.
104 '''
105 x1, x2 = self, self.others(other)
106 r1, r2, x = _rads3(rad1, rad2, radius)
107 if x:
108 x1, x2 = x2, x1
109 try:
110 n, q = x1.cross(x2), x1.dot(x2)
111 n2, q1 = n.length2, (_1_0 - q**2)
112 if n2 < EPS or isnear0(q1):
113 raise ValueError(_near_(_concentric_))
114 c1, c2 = cos(r1), cos(r2)
115 x0 = x1.times((c1 - q * c2) / q1).plus(
116 x2.times((c2 - q * c1) / q1))
117 n1 = _1_0 - x0.length2
118 if n1 < EPS:
119 raise ValueError(_too_(_distant_))
120 except ValueError as x:
121 raise IntersectionError(center=self, rad1=rad1,
122 other=other, rad2=rad2, cause=x)
123 n = n.times(sqrt(n1 / n2))
124 if n.length > EPS:
125 x1 = x0.plus(n)
126 x2 = x0.minus(n)
127 else: # abutting circles
128 x1 = x2 = x0
130 return (_xattrs(x1, self, _datum_, _name_),
131 _xattrs(x2, self, _datum_, _name_))
134class LatLonSphericalBase(LatLonBase):
135 '''(INTERNAL) Base class for spherical C{LatLon}s.
136 '''
137 _datum = Datums.Sphere # spherical L{Datum}
139 def __init__(self, latlonh, lon=None, height=0, datum=None, wrap=False, name=NN):
140 '''Create a spherical C{LatLon} point frome the given lat-, longitude and
141 height on the given datum.
143 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or
144 a previous C{LatLon} instance provided C{B{lon}=None}.
145 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or
146 C(None), indicating B{C{latlonh}} is a C{LatLon}.
147 @kwarg height: Optional height above (or below) the earth surface (C{meter},
148 same units as the datum's ellipsoid axes or radius).
149 @kwarg datum: Optional, spherical datum to use (L{Datum}, L{Ellipsoid},
150 L{Ellipsoid2}, L{a_f2Tuple}) or earth radius in C{meter},
151 conventionally).
152 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}
153 (C{bool}).
154 @kwarg name: Optional name (string).
156 @raise TypeError: If B{C{latlonh}} is not a C{LatLon} or B{C{datum}} not
157 spherical.
158 '''
159 LatLonBase.__init__(self, latlonh, lon=lon, height=height, wrap=wrap, name=name)
160 if datum not in (None, self.datum):
161 self.datum = datum
163 def bearingTo2(self, other, wrap=False, raiser=False):
164 '''Return the initial and final bearing (forward and reverse
165 azimuth) from this to an other point.
167 @arg other: The other point (C{LatLon}).
168 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
169 B{C{other}} point (C{bool}).
171 @return: A L{Bearing2Tuple}C{(initial, final)}.
173 @raise TypeError: The B{C{other}} point is not spherical.
175 @see: Methods C{initialBearingTo} and C{finalBearingTo}.
176 '''
177 # .initialBearingTo is inside .-Nvector and .-Trigonometry
178 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo
179 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser)
180 return Bearing2Tuple(i, f, name=self.name)
182 @property_doc_(''' this point's datum (L{Datum}).''')
183 def datum(self):
184 '''Get this point's datum (L{Datum}).
185 '''
186 return self._datum
188 @datum.setter # PYCHOK setter!
189 def datum(self, datum):
190 '''Set this point's datum I{without conversion} (L{Datum}, L{Ellipsoid},
191 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} spherical earth radius).
193 @raise TypeError: If B{C{datum}} invalid or not not spherical.
194 '''
195 d = _spherical_datum(datum, name=self.name, raiser=_datum_)
196 if self._datum != d:
197 _update_all(self)
198 self._datum = d
200 def finalBearingTo(self, other, wrap=False, raiser=False):
201 '''Return the final bearing (reverse azimuth) from this to
202 an other point.
204 @arg other: The other point (spherical C{LatLon}).
205 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
206 the B{C{other}} point (C{bool}).
208 @return: Final bearing (compass C{degrees360}).
210 @raise TypeError: The B{C{other}} point is not spherical.
212 @example:
214 >>> p = LatLon(52.205, 0.119)
215 >>> q = LatLon(48.857, 2.351)
216 >>> b = p.finalBearingTo(q) # 157.9
217 '''
218 p = self.others(other)
219 if wrap:
220 p = _unrollon(self, p, wrap=wrap)
221 # final bearing is the reverse of the other, initial one;
222 # .initialBearingTo is inside .-Nvector and .-Trigonometry
223 b = p.initialBearingTo(self, wrap=False, raiser=raiser)
224 return _umod_360(b + _180_0)
226 def intersecant2(self, circle, point, bearing, radius=R_M, exact=False,
227 height=None, wrap=False): # was=True
228 '''Compute the intersections of a circle and a line.
230 @arg circle: Radius of the circle centered at this location
231 (C{meter}, same units as B{C{radius}}) or a point
232 on the circle (this C{LatLon}).
233 @arg point: An other point in- or outside the circle (this C{LatLon}).
234 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360})
235 or a second point on the line (this C{LatLon}).
236 @kwarg radius: Mean earth radius (C{meter}, conventionally).
237 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
238 destination and distance, if C{False} use the basic
239 rhumb methods (C{bool}) or if C{None} use the I{great
240 circle} methods.
241 @kwarg height: Optional height for the intersection points (C{meter},
242 conventionally) or C{None}.
243 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
244 B{C{point}}, B{C{circle}} and/or B{C{bearing}} (C{bool}).
246 @return: 2-Tuple of the intersection points (representing a chord),
247 each an instance of this class. For a tangent line, each
248 point C{is} this very instance.
250 @raise IntersectionError: The circle and line do not intersect.
252 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}
253 or B{C{bearing}} invalid.
255 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
256 B{C{exact}} or B{C{height}}.
257 '''
258 p = self.others(point=point)
259 try:
260 return _intersecant2(self, circle, p, bearing, radius=radius, exact=exact,
261 height=height, wrap=wrap)
262 except (TypeError, ValueError) as x:
263 raise _xError(x, center=self, circle=circle, point=point, bearing=bearing,
264 exact=exact, wrap=wrap)
266 def maxLat(self, bearing):
267 '''Return the maximum latitude reached when travelling
268 on a great circle on given bearing from this point
269 based on Clairaut's formula.
271 The maximum latitude is independent of longitude
272 and the same for all points on a given latitude.
274 Negate the result for the minimum latitude (on the
275 Southern hemisphere).
277 @arg bearing: Initial bearing (compass C{degrees360}).
279 @return: Maximum latitude (C{degrees90}).
281 @raise ValueError: Invalid B{C{bearing}}.
282 '''
283 m = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi)))
284 return degrees90(m)
286 def minLat(self, bearing):
287 '''Return the minimum latitude reached when travelling
288 on a great circle on given bearing from this point.
290 @arg bearing: Initial bearing (compass C{degrees360}).
292 @return: Minimum latitude (C{degrees90}).
294 @see: Method L{maxLat} for more details.
296 @raise ValueError: Invalid B{C{bearing}}.
297 '''
298 return -self.maxLat(bearing)
300 def parse(self, strllh, height=0, sep=_COMMA_, name=NN):
301 '''Parse a string representing a similar, spherical C{LatLon}
302 point, consisting of C{"lat, lon[, height]"}.
304 @arg strllh: Lat, lon and optional height (C{str}),
305 see function L{pygeodesy.parse3llh}.
306 @kwarg height: Optional, default height (C{meter}).
307 @kwarg sep: Optional separator (C{str}).
308 @kwarg name: Optional instance name (C{str}),
309 overriding this name.
311 @return: The similar point (spherical C{LatLon}).
313 @raise ParseError: Invalid B{C{strllh}}.
314 '''
315 t = _MODS.dms.parse3llh(strllh, height=height, sep=sep)
316 r = self.classof(*t)
317 if name:
318 r.rename(name)
319 return r
321 @property_RO
322 def _radius(self):
323 '''(INTERNAL) Get this sphere's radius.
324 '''
325 return self.datum.ellipsoid.equatoradius
327 def _rhumbs3(self, other, wrap, r=False): # != .latlonBase._rhumbx3
328 '''(INTERNAL) Rhumb_ helper function.
330 @arg other: The other point (spherical C{LatLon}).
331 '''
332 p = self.others(other, up=2)
333 if wrap:
334 p = _unrollon(self, p, wrap=wrap)
335 a2, b2 = p.philam
336 a1, b1 = self.philam
337 # if |db| > 180 take shorter rhumb
338 # line across the anti-meridian
339 db = wrapPI(b2 - b1)
340 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1))
341 da = a2 - a1
342 if r:
343 # on Mercator projection, longitude distances shrink
344 # by latitude; the 'stretch factor' q becomes ill-
345 # conditioned along E-W line (0/0); use an empirical
346 # tolerance to avoid it
347 q = (da / dp) if fabs(dp) > EPS else cos(self.phi)
348 da = hypot(da, q * db) # angular distance radians
349 return da, db, dp
351 def rhumbAzimuthTo(self, other, radius=R_M, exact=False, wrap=False):
352 '''Return the azimuth (bearing) of a rhumb line (loxodrome)
353 between this and an other (spherical) point.
355 @arg other: The other point (spherical C{LatLon}).
356 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
357 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).
358 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
359 default C{False} for backward compatibility.
360 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
361 B{C{other}} point (C{bool}).
363 @return: Rhumb line azimuth (compass C{degrees180}).
365 @raise TypeError: The B{C{other}} point is incompatible or
366 B{C{radius}} is invalid.
368 @example:
370 >>> p = LatLon(51.127, 1.338)
371 >>> q = LatLon(50.964, 1.853)
372 >>> b = p.rhumbBearingTo(q) # 116.7
373 '''
374 if exact: # use series, always
375 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False,
376 radius=radius, wrap=wrap)
377 else:
378 _, db, dp = self._rhumbs3(other, wrap)
379 z = atan2d(db, dp) # see .rhumbx.Rhumb.Inverse
380 return z
382 @deprecated_method
383 def rhumbBearingTo(self, other): # unwrapped
384 '''DEPRECATED, use method C{.rhumbAzimuthTo}.'''
385 return wrap360(self.rhumbAzimuthTo(other)) # [0..360)
387 def rhumbDestination(self, distance, bearing, radius=R_M, height=None, exact=False):
388 '''Return the destination point having travelled the given distance
389 from this point along a rhumb line (loxodrome) at the given bearing.
391 @arg distance: Distance travelled (C{meter}, same units as B{C{radius}}),
392 may be negative if C{B{exact}=True}.
393 @arg bearing: Bearing (azimuth) at this point (compass C{degrees360}).
394 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
395 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if
396 C{B{exact}=True}.
397 @kwarg height: Optional height, overriding the default height
398 (C{meter}, same unit as B{C{radius}}).
399 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
400 default C{False} for backward compatibility.
402 @return: The destination point (spherical C{LatLon}).
404 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
405 B{C{radius}} or B{C{height}}.
407 @example:
409 >>> p = LatLon(51.127, 1.338)
410 >>> q = p.rhumbDestination(40300, 116.7) # 50.9642°N, 001.8530°E
411 '''
412 if exact: # use series, always
413 r = LatLonBase.rhumbDestination(self, distance, bearing, exact=False,
414 radius=radius, height=height)
415 else: # radius=None from .rhumbMidpointTo
416 if radius in (None, self._radius):
417 d, r = self.datum, radius
418 else:
419 d = _spherical_datum(radius, raiser=_radius_) # spherical only
420 r = d.ellipsoid.equatoradius
421 r = _angular(distance, r, low=_0_0) # distance=0 from .rhumbMidpointTo
423 a1, b1 = self.philam
424 sb, cb = sincos2(Bearing_(bearing))
426 da = r * cb
427 a2 = a1 + da
428 # normalize latitude if past pole
429 if a2 > PI_2:
430 a2 = PI - a2
431 elif a2 < -PI_2:
432 a2 = -PI - a2
434 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1))
435 # q becomes ill-conditioned on E-W course 0/0
436 q = (da / dp) if fabs(dp) > EPS else cos(a1)
437 b2 = (b1 + r * sb / q) if fabs(q) > EPS else b1
439 h = self._heigHt(height)
440 r = self.classof(degrees90(a2), degrees180(b2), datum=d, height=h)
441 return r
443 def rhumbDistanceTo(self, other, radius=R_M, exact=False, wrap=False):
444 '''Return the distance from this to an other point along
445 a rhumb line (loxodrome).
447 @arg other: The other point (spherical C{LatLon}).
448 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
449 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if
450 C{B{exact}=True}.
451 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
452 default C{False} for backward compatibility.
453 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
454 B{C{other}} point (C{bool}).
456 @return: Distance (C{meter}, the same units as B{C{radius}}
457 or C{radians} if B{C{radius}} is C{None}).
459 @raise TypeError: The B{C{other}} point is incompatible.
461 @raise ValueError: Invalid B{C{radius}}.
463 @example:
465 >>> p = LatLon(51.127, 1.338)
466 >>> q = LatLon(50.964, 1.853)
467 >>> d = p.rhumbDistanceTo(q) # 403100
468 '''
469 if exact: # use series, always
470 r = LatLonBase.rhumbDistanceTo(self, other, exact=False,
471 radius=radius, wrap=wrap)
472 if radius is None: # angular distance in radians
473 r = r / self._radius # /= chokes PyChecker
474 else:
475 # see <https://www.EdWilliams.org/avform.htm#Rhumb>
476 r, _, _ = self._rhumbs3(other, wrap, r=True)
477 if radius is not None:
478 r *= Radius(radius)
479 return r
481 def rhumbMidpointTo(self, other, height=None, radius=R_M, exact=False,
482 fraction=_0_5, wrap=False):
483 '''Return the (loxodromic) midpoint on the rhumb line between
484 this and an other point.
486 @arg other: The other point (spherical LatLon).
487 @kwarg height: Optional height, overriding the mean height
488 (C{meter}).
489 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
490 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).
491 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
492 default C{False} for backward compatibility.
493 @kwarg fraction: Midpoint location from this point (C{scalar}),
494 may be negative if C{B{exact}=True}.
495 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
496 B{C{other}} point (C{bool}).
498 @return: The (mid)point at the given B{C{fraction}} along
499 the rhumb line (spherical C{LatLon}).
501 @raise TypeError: The B{C{other}} point is incompatible.
503 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}
505 @example:
507 >>> p = LatLon(51.127, 1.338)
508 >>> q = LatLon(50.964, 1.853)
509 >>> m = p.rhumb_midpointTo(q)
510 >>> m.toStr() # '51.0455°N, 001.5957°E'
511 '''
512 if exact: # use series, always
513 r = LatLonBase.rhumbMidpointTo(self, other, exact=False,
514 radius=radius, height=height,
515 fraction=fraction, wrap=wrap)
516 elif fraction is not _0_5:
517 f = Scalar_(fraction=fraction) # low=_0_0
518 r, db, dp = self._rhumbs3(other, wrap, r=True) # radians
519 z = atan2b(db, dp)
520 h = self._havg(other, f=f, h=height)
521 r = self.rhumbDestination(r * f, z, radius=None, height=h)
523 else: # for backward compatibility, unwrapped
524 # see <https://MathForum.org/library/drmath/view/51822.html>
525 a1, b1 = self.philam
526 a2, b2 = self.others(other).philam
528 if fabs(b2 - b1) > PI:
529 b1 += PI2 # crossing anti-meridian
531 a3 = favg(a1, a2)
532 b3 = favg(b1, b2)
534 f1 = tanPI_2_2(a1)
535 if isnon0(f1):
536 f2 = tanPI_2_2(a2)
537 f = f2 / f1
538 if isnon0(f):
539 f = log(f)
540 if isnon0(f):
541 f3 = tanPI_2_2(a3)
542 b3 = fdot(map1(log, f1, f2, f3),
543 -b2, b1, b2 - b1) / f
545 d = self.datum if radius in (None, self._radius) else \
546 _spherical_datum(radius, name=self.name, raiser=_radius_)
547 h = self._havg(other, h=height)
548 r = self.classof(degrees90(a3), degrees180(b3), datum=d, height=h)
549 return r
551 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature
552 '''Convert this point to C{Nvector} components, I{including
553 height}.
555 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}
556 keyword arguments, ignored if
557 C{B{Nvector} is None}.
559 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}
560 if B{C{Nvector}} is C{None}.
562 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.
563 '''
564 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)
566 def toWm(self, radius=R_MA):
567 '''Convert this point to a I{WM} coordinate.
569 @kwarg radius: Optional earth radius (C{meter}).
571 @return: The WM coordinate (L{Wm}).
573 @see: Function L{pygeodesy.toWm} in module L{webmercator} for details.
574 '''
575 return _MODS.webmercator.toWm(self, radius=radius)
578def _intersecant2(c, r, p, b, radius=R_M, exact=False,
579 height=None, wrap=False):
580 # (INTERNAL) Intersect a circle and bearing, see L{intersecant2}
581 # above, separated to allow callers to embellish any exceptions
583 if wrap:
584 p = _unrollon(c, p, wrap=wrap)
585 nonexact = exact is None
587 if not isinstanceof(r, c.__class__, p.__class__):
588 r = Radius_(circle=r)
589 elif nonexact:
590 r = c.distanceTo(r, radius=radius, wrap=wrap)
591 elif isbool(exact):
592 r = c.rhumbDistanceTo(r, radius=radius, exact=exact, wrap=wrap)
593 else:
594 raise _ValueError(exact=exact)
596 if not isinstanceof(b, c.__class__, p.__class__):
597 b = Bearing(b)
598 elif nonexact:
599 b = p.initialBearingTo(b, wrap=wrap)
600 else:
601 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact, wrap=wrap)
603 d = p.distanceTo(c, radius=radius) if nonexact else \
604 p.rhumbDistanceTo(c, radius=radius, exact=exact)
605 if d > EPS:
606 a = p.initialBearingTo(c) if nonexact else wrap360(
607 p.rhumbAzimuthTo(c, radius=radius, exact=exact))
608 s, c = sincos2d(b - a)
609 s = sqrt_a(r, fabs(s * d))
610 if s > r:
611 raise IntersectionError(_too_(Fmt.distant(s)))
612 elif (r - s) < EPS:
613 return p, p # tangent
614 c *= d
615 else: # coincindent
616 s, c = r, 0
618 a = b + _180_0
619 if nonexact:
620 b = p.destination(s + c, b, radius=radius, height=height)
621 a = p.destination(s - c, a, radius=radius, height=height)
622 else:
623 b = p.rhumbDestination(s + c, b, radius=radius, height=height, exact=exact)
624 a = p.rhumbDestination(s - c, a, radius=radius, height=height, exact=exact)
625 return b, a # in bearing direction first
628__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase)
630# **) MIT License
631#
632# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
633#
634# Permission is hereby granted, free of charge, to any person obtaining a
635# copy of this software and associated documentation files (the "Software"),
636# to deal in the Software without restriction, including without limitation
637# the rights to use, copy, modify, merge, publish, distribute, sublicense,
638# and/or sell copies of the Software, and to permit persons to whom the
639# Software is furnished to do so, subject to the following conditions:
640#
641# The above copyright notice and this permission notice shall be included
642# in all copies or substantial portions of the Software.
643#
644# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
645# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
646# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
647# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
648# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
649# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
650# OTHER DEALINGS IN THE SOFTWARE.