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