1
2# -*- coding: utf-8 -*-
3
4u'''(INTERNAL) Ellipsoidal direct/inverse geodesy base class
5C{LatLonEllipsoidalBaseDI} and functions.
6'''
7# make sure int/int division yields float quotient, see .basics
8from __future__ import division
9
10from pygeodesy.basics import isscalar, issubclassof
11from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, Property_RO
12from pygeodesy.errors import _AssertionError, IntersectionError, _IsnotError, \
13                             _ValueError, _xellipsoidal, _xError, _xkwds_not
14from pygeodesy.fmath import favg, fmean_, fsum_
15from pygeodesy.formy import opposing, _radical2
16from pygeodesy.interns import EPS, PI, PI_4, _concentric_, _datum_, _epoch_, \
17                             _exceed_PI_radians_, _height_, _no_, _near_, \
18                             _reframe_, _too_, _0_0, _1_5, _3_0
19from pygeodesy.lazily import _ALL_DOCS
20from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \
21                           Intersection3Tuple, _LL4Tuple
22# from pygeodesy.props import Property_RO  # from .ellipsoidalBase
23from pygeodesy.streprs import Fmt
24from pygeodesy.units import Height, Radius_, Scalar, _1mm as _TOL_M
25from pygeodesy.utily import m2degrees, unroll180, wrap90, wrap180, wrap360
26
27__all__ = ()
28__version__ = '21.09.16'
29
30_B2END = _1_5  # _intersect3 bearing to end point
31_TRIPS = 17    # _intersect3, _intersects2, _nearestOn interations, 6 is sufficient
32
33
34class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase):
35    '''(INTERNAL) Base class for C{ellipsoidal*.LatLon} classes
36       with I{overloaded} C{Direct} and C{Inverse} methods.
37    '''
38
39    def bearingTo2(self, other, wrap=False):
40        '''Compute the initial and final bearing (forward and reverse
41           azimuth) from this to an other point, using this C{Inverse}
42           method.  See methods L{initialBearingTo} and L{finalBearingTo}
43           for more details.
44
45           @arg other: The other point (C{LatLon}).
46           @kwarg wrap: Wrap and unroll longitudes (C{bool}).
47
48           @return: A L{Bearing2Tuple}C{(initial, final)}.
49
50           @raise TypeError: The B{C{other}} point is not C{LatLon}.
51
52           @raise ValueError: If this and the B{C{other}} point's L{Datum}
53                              ellipsoids are not compatible.
54        '''
55        r = self._Inverse(other, wrap)
56        return Bearing2Tuple(r.initial, r.final, name=self.name)
57
58    def destination(self, distance, bearing, height=None):
59        '''Compute the destination point after having travelled for
60           the given distance from this point along a geodesic given
61           by an initial bearing, using this C{Direct} method.  See
62           method L{destination2} for more details.
63
64           @arg distance: Distance (C{meter}).
65           @arg bearing: Initial bearing in (compass C{degrees360}).
66           @kwarg height: Optional height, overriding the default
67                          height (C{meter}, same units as C{distance}).
68
69           @return: The destination point (C{LatLon}).
70        '''
71        return self._Direct(distance, bearing, self.classof, height).destination
72
73    def destination2(self, distance, bearing, height=None):
74        '''Compute the destination point and the final bearing (reverse
75           azimuth) after having travelled for the given distance from
76           this point along a geodesic given by an initial bearing,
77           using this C{Direct} method.
78
79           The distance must be in the same units as this point's datum
80           axes, conventionally C{meter}.  The distance is measured on
81           the surface of the ellipsoid, ignoring this point's height.
82
83           The initial and final bearing (forward and reverse azimuth)
84           are in compass C{degrees360}.
85
86           The destination point's height and datum are set to this
87           point's height and datum, unless the former is overridden.
88
89           @arg distance: Distance (C{meter}).
90           @arg bearing: Initial bearing (compass C{degrees360}).
91           @kwarg height: Optional height, overriding the default
92                          height (C{meter}, same units as C{distance}).
93
94           @return: A L{Destination2Tuple}C{(destination, final)}.
95        '''
96        r = self._Direct(distance, bearing, self.classof, height)
97        return self._xnamed(r)
98
99    def _Direct(self, distance, bearing, LL, height):  # overloaded by I{Vincenty}
100        '''(INTERNAL) I{Karney}'s C{Direct} method.
101
102           @return: A L{Destination2Tuple}C{(destination, final)} or
103                    a L{Destination3Tuple}C{(lat, lon, final)} if
104                    B{C{LL}} is C{None}.
105        '''
106        g = self.geodesic
107        r = g.Direct3(self.lat, self.lon, bearing, distance)
108        if LL:
109            r = self._Direct2Tuple(LL, height, r)
110        return r
111
112    def _Direct2Tuple(self, LL, height, r):
113        '''(INTERNAL) Helper for C{._Direct} result L{Destination2Tuple}.
114        '''
115        h = self.height if height is None else height
116        d = LL(wrap90(r.lat), wrap180(r.lon), height=h, datum=self.datum, name=self.name,
117                 **_xkwds_not(None, epoch=self.epoch, reframe=self.reframe))
118        return Destination2Tuple(d, wrap360(r.final))
119
120    def distanceTo(self, other, wrap=False, **unused):  # ignore radius=R_M
121        '''Compute the distance between this and an other point
122           along a geodesic, using this C{Inverse} method. See method
123           L{distanceTo3} for more details.
124
125           @arg other: The other point (C{LatLon}).
126           @kwarg wrap: Wrap and unroll longitudes (C{bool}).
127
128           @return: Distance (C{meter}).
129
130           @raise TypeError: The B{C{other}} point is not C{LatLon}.
131
132           @raise ValueError: If this and the B{C{other}} point's L{Datum}
133                              ellipsoids are not compatible.
134        '''
135        return self._Inverse(other, wrap, azis=False).distance
136
137    def distanceTo3(self, other, wrap=False):
138        '''Compute the distance, the initial and final bearing along
139           a geodesic between this and an other point, using this
140           C{Inverse} method.
141
142           The distance is in the same units as this point's datum axes,
143           conventionally meter.  The distance is measured on the surface
144           of the ellipsoid, ignoring this point's height.
145
146           The initial and final bearing (forward and reverse azimuth)
147           are in compass C{degrees360} from North.
148
149           @arg other: Destination point (C{LatLon}).
150           @kwarg wrap: Wrap and unroll longitudes (C{bool}).
151
152           @return: A L{Distance3Tuple}C{(distance, initial, final)}.
153
154           @raise TypeError: The B{C{other}} point is not C{LatLon}.
155
156           @raise ValueError: If this and the B{C{other}} point's L{Datum}
157                              ellipsoids are not compatible.
158        '''
159        return self._xnamed(self._Inverse(other, wrap))
160
161    def finalBearingOn(self, distance, bearing):
162        '''Compute the final bearing (reverse azimuth) after having
163           travelled for the given distance along a geodesic given by
164           an initial bearing from this point, using this C{Direct}
165           method.  See method L{destination2} for more details.
166
167           @arg distance: Distance (C{meter}).
168           @arg bearing: Initial bearing (compass C{degrees360}).
169
170           @return: Final bearing (compass C{degrees360}).
171        '''
172        return self._Direct(distance, bearing, None, None).final
173
174    def finalBearingTo(self, other, wrap=False):
175        '''Compute the final bearing (reverse azimuth) after having
176           travelled along a geodesic from this point to an other
177           point, using this C{Inverse} method.  See method
178           L{distanceTo3} for more details.
179
180           @arg other: The other point (C{LatLon}).
181           @kwarg wrap: Wrap and unroll longitudes (C{bool}).
182
183           @return: Final bearing (compass C{degrees360}).
184
185           @raise TypeError: The B{C{other}} point is not C{LatLon}.
186
187           @raise ValueError: If this and the B{C{other}} point's L{Datum}
188                              ellipsoids are not compatible.
189        '''
190        return self._Inverse(other, wrap).final
191
192    @Property_RO
193    def geodesic(self):  # overloaded by I{Karney}'s, N/A for I{Vincenty}
194        '''N/A, invalid (C{None} I{always}).
195        '''
196        return None  # PYCHOK no cover
197
198    def initialBearingTo(self, other, wrap=False):
199        '''Compute the initial bearing (forward azimuth) to travel
200           along a geodesic from this point to an other point,
201           using this C{Inverse} method.  See method L{distanceTo3}
202           for more details.
203
204           @arg other: The other point (C{LatLon}).
205           @kwarg wrap: Wrap and unroll longitudes (C{bool}).
206
207           @return: Initial bearing (compass C{degrees360}).
208
209           @raise TypeError: The B{C{other}} point is not C{LatLon}.
210
211           @raise ValueError: If this and the B{C{other}} point's L{Datum}
212                              ellipsoids are not compatible.
213        '''
214        return self._Inverse(other, wrap).initial
215
216    def intermediateTo(self, other, fraction, height=None, wrap=False):
217        '''Return the point at given fraction along the geodesic between
218           this and an other point, using this C{Direct} and C{Inverse}
219           methods.
220
221           @arg other: The other point (C{LatLon}).
222           @arg fraction: Fraction between both points ranging from
223                          0, meaning this to 1, the other point (C{float}).
224           @kwarg height: Optional height, overriding the fractional
225                          height (C{meter}).
226           @kwarg wrap: Wrap and unroll longitudes (C{bool}).
227
228           @return: Intermediate point (C{LatLon}).
229
230           @raise TypeError: The B{C{other}} point is not C{LatLon}.
231
232           @raise UnitError: Invalid B{C{fraction}} or B{C{height}}.
233
234           @raise ValueError: If this and the B{C{other}} point's L{Datum}
235                              ellipsoids are not compatible.
236
237           @see: Methods L{distanceTo3} and L{destination}.
238        '''
239        t = self.distanceTo3(other, wrap=wrap)
240        f = Scalar(fraction=fraction)
241        h = self._havg(other, f=f) if height is None else Height(height)
242        return self.destination(t.distance * f, t.initial, height=h)
243
244    def _Inverse(self, other, wrap, **unused):  # azis=False, overloaded by I{Vincenty}
245        '''(INTERNAL) I{Karney}'s C{Inverse} method.
246
247           @return: A L{Distance3Tuple}C{(distance, initial, final)}.
248
249           @raise TypeError: The B{C{other}} point is not L{LatLon}.
250
251           @raise ValueError: If this and the B{C{other}} point's
252                              L{Datum} ellipsoids are not compatible.
253        '''
254        _ = self.ellipsoids(other)
255        g = self.geodesic
256        _, lon = unroll180(self.lon, other.lon, wrap=wrap)
257        return g.Inverse3(self.lat, self.lon, other.lat, lon)
258
259
260def _Equidistant00(equidistant, p1):
261    '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.
262    '''
263    if equidistant is None or not callable(equidistant):
264        equidistant = p1.Equidistant
265    else:
266        from pygeodesy.azimuthal import _Equidistants
267        if not issubclassof(equidistant, *_Equidistants):  # PYCHOK no cover
268            t = tuple(_.__name__ for _ in _Equidistants)
269            raise _IsnotError(*t, equidistant=equidistant)
270    return equidistant(0, 0, p1.datum)
271
272
273def _intersect3(s1, end1, s2, end2, height=None, wrap=True,  # MCCABE 16
274                equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
275    '''(INTERNAL) Intersect two (ellipsoidal) path, see ellipsoidal method
276       L{intersection3}, separated to allow callers to embellish any exceptions.
277    '''
278    from pygeodesy.sphericalTrigonometry import _intersect as _si, LatLon as _LLS
279    from pygeodesy.vector3d import _intersect3d3 as _vi3
280
281    def _b_d(s, e, w, t, h=_0_0):
282        # compute opposing and distance
283        t = s.classof(t.lat, t.lon, height=h, name=t.name)
284        t = s.distanceTo3(t, wrap=w)  # Distance3Tuple
285        b = opposing(e, t.initial)  # "before" start
286        return b, t.distance
287
288    def _b_e(s, e, w, t):
289        # compute an end point along the initial bearing
290        # about 1.5 times the distance to the gu-/estimate, at
291        # least 1/8 and at most 3/8 of the earth perimeter like
292        # radians in .sphericalTrigonometry._int3d2 and bearing
293        # comparison in .sphericalTrigonometry._intb
294        b, d = _b_d(s, e, w, t, h=t.height)
295        m = s.ellipsoid().R2x * PI_4  # authalic exact
296        d = min(max(d * _B2END, m), m * _3_0)
297        e = s.destination(d, e)
298        return b, (_unrollon(s, e) if w else e)
299
300    def _e_ll(s, e, w, **end):
301        # return 2-tuple (end, False if bearing else True)
302        ll = not isscalar(e)
303        if ll:
304            e = s1.others(**end)
305            if w:  # unroll180 == .karney._unroll2
306                e = _unrollon(s, e)
307        return e, ll
308
309    def _o(o, b, n, s, t, eps):
310        # determine C{o}utside before, on or after start point
311        if not o:  # intersection may be on start
312            from pygeodesy.latlonBase import _isequalTo
313            if _isequalTo(s, t, eps=eps):
314                return o
315        return -n if b else n
316
317    # E = s1.ellipsoids(s2)
318    # assert E == s1.ellispoids(e1) == s1.ellipsoids(e2)
319
320    e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)
321    e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)
322
323    # get the azimuthal equidistant projection
324    A = _Equidistant00(equidistant, s1)
325    e = max(tol, EPS)
326
327    # gu-/estimate initial intersection, spherically ...
328    t = _si(_LLS(s1.lat, s1.lon, height=s1.height),
329           (_LLS(e1.lat, e1.lon, height=e1.height) if ll1 else e1),
330            _LLS(s2.lat, s2.lon, height=s2.height),
331           (_LLS(e2.lat, e2.lon, height=e2.height) if ll2 else e2),
332            height=height, wrap=False, LatLon=_LLS)  # unrolled already
333    h, n = t.height, t.name
334
335    if not ll1:
336        b1, e1 = _b_e(s1, e1, wrap, t)
337    if not ll2:
338        b2, e2 = _b_e(s2, e2, wrap, t)
339
340    # ... and iterate as Karney describes, @see:
341    # LatLonEllipsoidalBase.LatLon.intersections2
342    c = m = None  # force first d == c to False
343    for i in range(1, _TRIPS):
344        A.reset(t.lat, t.lon)  # gu-/estimate as origin
345        # convert start and end points to projection
346        # space and compute an intersection there
347        v, o1, o2 = _vi3(A.forward(s1.lat, s1.lon),
348                         A.forward(e1.lat, e1.lon),
349                         A.forward(s2.lat, s2.lon),
350                         A.forward(e2.lat, e2.lon),
351                         eps=e, useZ=False)
352        # convert intersection back to geodetic
353        t, d = A._reverse2(v.x, v.y)
354        # break if below tolerance or if unchanged
355        if d < e or d == c:
356            t._iteration = i
357            break
358        if m is None or m > d:
359            m = d
360        c = d
361    else:
362        raise IntersectionError(_no_(Fmt.convergence(m)))
363
364    # like .sphericalTrigonometry._intersect, if this intersection
365    # is "before" the first point, use the antipodal intersection
366    if not (ll1 or ll2):  # end1 and end2 are an initial bearing
367        b1, _ = _b_d(s1, end1, wrap, t)
368        if b1:
369            t  = t.antipodal()
370            b1 = not b1
371        b2, _ = _b_d(s2, end2, wrap, t)
372
373    r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1, name=n)
374    r._iteration = t._iteration  # _NamedTuple._iteration
375    return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)),
376                                 (o2 if ll2 else _o(o2, b2, 2, s2, t, e)))
377
378
379def _intersection3(start1, end1, start2, end2, height=None, wrap=True,
380                   equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
381    '''(INTERNAL) Iteratively compute the intersection point of two paths,
382       each defined by two (ellipsoidal) points or an (ellipsoidal) start
383       point and an initial bearing from North.
384    '''
385    s1 = _xellipsoidal(start1=start1)
386    s2 =  s1.others(start2=start2)
387    try:
388        return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap,
389                                          equidistant=equidistant, tol=tol,
390                                               LatLon=LatLon, **LatLon_kwds)
391    except (TypeError, ValueError) as x:
392        raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
393
394
395def _intersections2(center1, radius1, center2, radius2, height=None, wrap=True,
396                    equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
397    '''(INTERNAL) Iteratively compute the intersection points of two circles,
398       each defined by an (ellipsoidal) center point and a radius.
399    '''
400    c1 = _xellipsoidal(center1=center1)
401    c2 = c1.others(center2=center2)
402    try:
403        return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap,
404                                                 equidistant=equidistant, tol=tol,
405                                                      LatLon=LatLon, **LatLon_kwds)
406    except (TypeError, ValueError) as x:
407        raise _xError(x, center1=center1, radius1=radius1,
408                         center2=center2, radius2=radius2)
409
410
411def _intersects2(c1, radius1, c2, radius2, height=None, wrap=True,  # MCCABE 16
412                 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
413    '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2}
414       above, separated to allow callers to embellish any exceptions.
415    '''
416    from pygeodesy.sphericalTrigonometry import _intersects2 as _si2, LatLon as _LLS
417    from pygeodesy.vector3d import _intersects2 as _vi2
418
419    def _latlon4(t, h, n, c):
420        r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c, name=n)
421        r._iteration = t.iteration  # ._iteration for tests
422        return r
423
424    r1 = Radius_(radius1=radius1)
425    r2 = Radius_(radius2=radius2)
426
427    E = c1.ellipsoids(c2)
428    # get the azimuthal equidistant projection
429    A = _Equidistant00(equidistant, c1)
430
431    if r1 < r2:
432        c1, c2 = c2, c1
433        r1, r2 = r2, r1
434
435    if r1 > (min(E.b, E.a) * PI):
436        raise _ValueError(_exceed_PI_radians_)
437
438    if wrap:  # unroll180 == .karney._unroll2
439        c2 = _unrollon(c1, c2)
440
441    # distance between centers and radii are
442    # measured along the ellipsoid's surface
443    m = c1.distanceTo(c2, wrap=False)  # meter
444    if m < max(r1 - r2, EPS):
445        raise IntersectionError(_near_(_concentric_))
446    if fsum_(r1, r2, -m) < 0:
447        raise IntersectionError(_too_(Fmt.distant(m)))
448
449    f = _radical2(m, r1, r2).ratio  # "radical fraction"
450    r = E.rocMean(favg(c1.lat, c2.lat, f=f))
451    e = max(m2degrees(tol, radius=r), EPS)
452
453    # gu-/estimate initial intersections, spherically ...
454    t1, t2 = _si2(_LLS(c1.lat, c1.lon, height=c1.height), r1,
455                  _LLS(c2.lat, c2.lon, height=c2.height), r2,
456                   radius=r, height=height, wrap=False, too_d=m)  # unrolled already
457    h, n = t1.height, t1.name
458
459    # ... and iterate as Karney describes, @see:
460    # LatLonEllipsoidalBase.LatLon.intersections2
461    ts, ta = [], None
462    for t in ((t1,) if t1 is t2 else (t1, t2)):
463        c = m = None  # force first d == c to False
464        for i in range(1, _TRIPS):
465            A.reset(t.lat, t.lon)  # gu-/estimate as origin
466            # convert centers to projection space
467            t1 = A.forward(c1.lat, c1.lon)
468            t2 = A.forward(c2.lat, c2.lon)
469            # compute intersections in projection space
470            v1, v2 = _vi2(t1, r1,  # XXX * t1.scale?,
471                          t2, r2,  # XXX * t2.scale?,
472                          sphere=False, too_d=m)
473            # convert intersections back to geodetic
474            t1, d1 = A._reverse2(v1.x, v1.y)
475            if v1 is v2:  # abutting
476                t, d = t1, d1  # PYCHOK no cover
477            else:
478                t2, d2 = A._reverse2(v2.x, v2.y)
479                # consider only the closer intersection
480                t, d = (t1, d1) if d1 < d2 else (t2, d2)
481            # break if below tolerance or if unchanged
482            if d < e or d == c:
483                t._iteration = i  # _NamedTuple._iteration
484                ts.append(t)
485                if v1 is v2:  # abutting
486                    ta = t  # PYCHOK no coves
487                break
488            c = d
489            if m is None or m > d:
490                m = d
491        else:
492            raise IntersectionError(_no_(Fmt.convergence(m)))
493
494    if ta:  # abutting circles
495        pass  # PYCHOK no cover
496    elif len(ts) == 2:
497        return (_latlon4(ts[0], h, n, c1),
498                _latlon4(ts[1], h, n, c2))
499    elif len(ts) == 1:  # PYCHOK no cover
500        ta = ts[0]  # assume abutting
501    else:  # PYCHOK no cover
502        raise _AssertionError(ts=ts)
503    r = _latlon4(ta, h, n, c1)
504    return r, r
505
506
507def _nearestOn(point, point1, point2, within=True, height=None, wrap=True,
508               equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
509    '''(INTERNAL) Get closest point, imported by .ellipsoidalExact,
510       -GeodSolve, -Karney and -Vincenty to embellish exceptions.
511    '''
512    try:
513        p = _xellipsoidal(point=point)
514        return _nearestOne(p, point1, point2, within=within, height=height, wrap=wrap,
515                           equidistant=equidistant, tol=tol, LatLon=LatLon, **LatLon_kwds)
516    except (TypeError, ValueError) as x:
517        raise _xError(x, point=point, point1=point1, point2=point2)
518
519
520def _nearestOne(p, point1, point2, within=True, height=None, wrap=True,
521                equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
522    '''(INTERNAL) Get closest point, like L{_intersects2} above,
523       separated to allow callers to embellish any exceptions.
524    '''
525    from pygeodesy.sphericalNvector import LatLon as _LLS
526    from pygeodesy.vector2d import _nearestOn as _vnOn, Vector3d
527
528    def _v(t, h):
529        return Vector3d(t.x, t.y, h)
530
531    p1 = p.others(point1=point1)
532    p2 = p.others(point2=point2)
533
534    _ = p.ellipsoids(p1)
535    E = p.ellipsoids(p2)
536
537    if wrap:
538        p1 = _unrollon(p,  p1)
539        p2 = _unrollon(p,  p2)
540        p2 = _unrollon(p1, p2)
541
542    r = E.rocMean(fmean_(p.lat, p1.lat, p2.lat))
543    e = max(m2degrees(tol, radius=r), EPS)
544
545    # get the azimuthal equidistant projection
546    A = _Equidistant00(equidistant, p)
547
548    # gu-/estimate initial nearestOn, spherically ... wrap=False, only!
549    t = _LLS(p.lat,  p.lon,  height=p.height).nearestOn(
550        _LLS(p1.lat, p1.lon, height=p1.height),
551        _LLS(p2.lat, p2.lon, height=p2.height), within=within, height=height)
552    n = t.name
553    if height is False:  # PYCHOK no cover
554        h  = t.height  # use heights as Z
555        h1 = p1.height
556        h2 = p2.height
557    else:
558        h = h1 = h2 = _0_0
559
560    # ... and iterate as Karney describes, @see:
561    # LatLonEllipsoidalBase.LatLon.intersections2
562    c = m = None  # force first d == c to False
563    # closest to origin, .z to interpolate height
564    p = Vector3d(0, 0, h)
565    for i in range(1, _TRIPS):
566        A.reset(t.lat, t.lon)  # gu-/estimate as origin
567        # convert points to projection space
568        # and compute the nearest one there
569        v = _vnOn(p, _v(A.forward(p1.lat, p1.lon), h1),
570                     _v(A.forward(p2.lat, p2.lon), h2),
571                      within=within)
572        # convert nearest one back to geodetic
573        t, d = A._reverse2(v.x, v.y)
574        # break if below tolerance or if unchanged
575        if d < e or d == c:
576            t._iteration = i  # _NamedTuple._iteration
577            if height is False:
578                h = v.z  # nearest interpolated
579            break
580        c = d
581        if m is None or m > d:
582            m = d
583    else:
584        raise _ValueError(_no_(Fmt.convergence(m)))
585
586    r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p, name=n)
587    r._iteration = t.iteration  # ._iteration for tests
588    return r
589
590
591def _unrollon(p1, p2):  # unroll180 == .karney._unroll2
592    # wrap, unroll and replace longitude if different
593    _, lon = unroll180(p1.lon, p2.lon, wrap=True)
594    if abs(lon - p2.lon) > EPS:
595        p2 = p2.classof(p2.lat, wrap180(lon), **_xkwds_not(None,
596                             height=getattr(p2, _height_,  None),
597                              datum=getattr(p2, _datum_,   None),
598                              epoch=getattr(p2, _epoch_,   None),
599                            reframe=getattr(p2, _reframe_, None)))  # PYCHOK indent
600    return p2
601
602
603__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI)
604
605# **) MIT License
606#
607# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved.
608#
609# Permission is hereby granted, free of charge, to any person obtaining a
610# copy of this software and associated documentation files (the "Software"),
611# to deal in the Software without restriction, including without limitation
612# the rights to use, copy, modify, merge, publish, distribute, sublicense,
613# and/or sell copies of the Software, and to permit persons to whom the
614# Software is furnished to do so, subject to the following conditions:
615#
616# The above copyright notice and this permission notice shall be included
617# in all copies or substantial portions of the Software.
618#
619# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
620# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
621# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
622# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
623# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
624# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
625# OTHER DEALINGS IN THE SOFTWARE.
626