
# -*- coding: utf-8 -*-

u'''Trigonometric classes geodetic (lat-/longitude) L{LatLon} and
geocentric (ECEF) L{Cartesian} and functionsL{areaOf}, L{intersection},
L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn2} and L{perimeterOf},
I{all spherical}.

Pure Python implementation of geodetic (lat-/longitude) methods using
spherical trigonometry, transcribed from JavaScript originals by
I{(C) Chris Veness 2011-2016} published under the same MIT Licence**, see
U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.

@newfield example: Example, Examples
'''

from pygeodesy.datum import R_M
from pygeodesy.fmath import EPS, acos1, favg, fdot, fmean, fsum, \
                            isscalar, map1
from pygeodesy.formy import antipode, bearing_, haversine_
from pygeodesy.lazily import _ALL_LAZY, _2kwds
from pygeodesy.named import LatLon3Tuple, NearestOn3Tuple, _xnamed
from pygeodesy.points import _imdex2, ispolar, nearestOn5 as _nearestOn5
from pygeodesy.sphericalBase import CartesianSphericalBase, LatLonSphericalBase
from pygeodesy.utily import PI2, PI_2, PI_4, degrees90, degrees180, \
                            degrees2m, iterNumpy2, radiansPI2, \
                            sincos2, tan_2, unrollPI, wrapPI
from pygeodesy.vector3d import CrossError, crosserrors, Vector3d, sumOf

from math import asin, atan2, copysign, cos, degrees, hypot, \
                 radians, sin

# all public contants, classes and functions
__all__ = _ALL_LAZY.sphericalTrigonometry + (
          'Cartesian', 'LatLon',  # classes
          'areaOf',  # functions
          'intersection', 'ispolar', 'isPoleEnclosedBy',  # DEPRECATED, use ispolar
          'meanOf',
          'nearestOn2', 'nearestOn3',
          'perimeterOf',
          'sumOf')  # == vector3d.sumOf
__version__ = '19.10.21'


class Cartesian(CartesianSphericalBase):
    '''Extended to convert geocentric, L{Cartesian} points to
       spherical, geodetic L{LatLon}.
    '''

    def toLatLon(self, **kwds):  # PYCHOK LatLon=LatLon
        '''Convert this cartesian point to an C{Nvector}-based
           geodetic point.

           @keyword kwds: Optional, additional B{C{LatLon}} keyword
                          arguments, ignored if C{B{LatLon}=None}.
                          For example, use C{LatLon=...} to override
                          the L{LatLon} (sub-)class or specify

           @return: The B{C{LatLon}} point (L{LatLon}) or if
                    C{B{LatLon}=None}, an L{Ecef9Tuple}C{(x, y, z,
                    lat, lon, height, C, M, datum)} with C{C} and
                    C{M} if available.

           @raise TypeError: Invalid B{C{LatLon}}, B{C{datum}}
                             or B{C{kwds}}.
        '''
        kwds = _2kwds(kwds, LatLon=LatLon, datum=self.datum)
        return CartesianSphericalBase.toLatLon(self, **kwds)


class LatLon(LatLonSphericalBase):
    '''New point on spherical model earth model.

       @example:

       >>> p = LatLon(52.205, 0.119)  # height=0
    '''

    def _trackDistanceTo3(self, start, end, radius, wrap):
        '''(INTERNAL) Helper for along-/crossTrackDistanceTo.
        '''
        self.others(start, name='start')
        self.others(end, name='end')

        r = float(radius)
        if r < EPS:
            raise ValueError('%s invalid: %r' % ('radius', radius))

        r = start.distanceTo(self, r, wrap=wrap) / r
        b = radians(start.initialBearingTo(self, wrap=wrap))
        e = radians(start.initialBearingTo(end, wrap=wrap))

        x = asin(sin(r) * sin(b - e))
        return r, x, e - b

    def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
        '''Compute the (signed) distance from the start to the closest
           point on the great circle path defined by a start and an
           end point.

           That is, if a perpendicular is drawn from this point to the
           great circle path, the along-track distance is the distance
           from the start point to the point where the perpendicular
           crosses the path.

           @param start: Start point of great circle path (L{LatLon}).
           @param end: End point of great circle path (L{LatLon}).
           @keyword radius: Optional, mean earth radius (C{meter}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: Distance along the great circle path (C{meter},
                    same units as B{C{radius}}), positive if after the
                    B{C{start}} toward the B{C{end}} point of the path or
                    negative if before the B{C{start}} point.

           @raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}.

           @example:

           >>> p = LatLon(53.2611, -0.7972)

           >>> s = LatLon(53.3206, -1.7297)
           >>> e = LatLon(53.1887, 0.1334)
           >>> d = p.alongTrackDistanceTo(s, e)  # 62331.58
        '''
        r, x, b = self._trackDistanceTo3(start, end, radius, wrap)
        cx = cos(x)
        if abs(cx) > EPS:
            return copysign(acos1(cos(r) / cx), cos(b)) * radius
        else:
            return 0.0

    def bearingTo(self, other, wrap=False, raiser=False):
        '''DEPRECATED, use method C{initialBearingTo}.
        '''
        return self.initialBearingTo(other, wrap=wrap, raiser=raiser)

    def crossingParallels(self, other, lat, wrap=False):
        '''Return the pair of meridians at which a great circle defined
           by this and an other point crosses the given latitude.

           @param other: The other point defining great circle (L{LatLon}).
           @param lat: Latitude at the crossing (C{degrees}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
                    C{None} if the great circle doesn't reach B{C{lat}}.
        '''
        self.others(other)

        a1, b1 = self.to2ab()
        a2, b2 = other.to2ab()

        a = radians(lat)
        db, b2 = unrollPI(b1, b2, wrap=wrap)

        sa,  ca,  sa1, ca1, \
        sa2, ca2, sdb, cdb = sincos2(a, a1, a2, db)

        x = sa1 * ca2 * ca * sdb
        y = sa1 * ca2 * ca * cdb - ca1 * sa2 * ca
        z = ca1 * ca2 * sa * sdb

        h = hypot(x, y)
        if h < EPS or abs(z) > h:
            return None  # great circle doesn't reach latitude

        m = atan2(-y, x) + b1  # longitude at max latitude
        d = acos1(z / h)  # delta longitude to intersections
        return degrees180(m - d), degrees180(m + d)

    def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
        '''Compute the (signed) distance from this point to the great
           circle defined by a start and an end point.

           @param start: Start point of great circle path (L{LatLon}).
           @param end: End point of great circle path (L{LatLon}).
           @keyword radius: Optional, mean earth radius (C{meter}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: Distance to great circle (negative if to the
                    left or positive if to the right of the path).

           @raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}.

           @example:

           >>> p = LatLon(53.2611, -0.7972)

           >>> s = LatLon(53.3206, -1.7297)
           >>> e = LatLon(53.1887, 0.1334)
           >>> d = p.crossTrackDistanceTo(s, e)  # -307.5
        '''
        _, x, _ = self._trackDistanceTo3(start, end, radius, wrap)
        return x * radius

    def destination(self, distance, bearing, radius=R_M, height=None):
        '''Locate the destination from this point after having
           travelled the given distance on the given initial bearing.

           @param distance: Distance travelled (C{meter}, same units as
                            B{C{radius}}).
           @param bearing: Bearing from this point (compass C{degrees360}).
           @keyword radius: Optional, mean earth radius (C{meter}).
           @keyword height: Optional height at destination (C{meter},
                            same units a B{C{radius}}).

           @return: Destination point (L{LatLon}).

           @example:

           >>> p1 = LatLon(51.4778, -0.0015)
           >>> p2 = p1.destination(7794, 300.7)
           >>> p2.toStr()  # '51.5135°N, 000.0983°W'

           @JSname: I{destinationPoint}.
        '''
        a, b = self.to2ab()

        r = float(distance) / float(radius)  # angular distance in radians
        t = radians(bearing)

        a, b = _destination2(a, b, r, t)
        h = self.height if height is None else height
        return self.classof(a, b, height=h)

    def distanceTo(self, other, radius=R_M, wrap=False):
        '''Compute the distance from this to an other point.

           @param other: The other point (L{LatLon}).
           @keyword radius: Optional, mean earth radius (C{meter}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: Distance between this and the B{C{other}} point
                    (C{meter}, same units as B{C{radius}}).

           @raise TypeError: The B{C{other}} point is not L{LatLon}.

           @example:

           >>> p1 = LatLon(52.205, 0.119)
           >>> p2 = LatLon(48.857, 2.351);
           >>> d = p1.distanceTo(p2)  # 404300
        '''
        self.others(other)

        a1, b1 = self.to2ab()
        a2, b2 = other.to2ab()

        db, b2 = unrollPI(b1, b2, wrap=wrap)
        r = haversine_(a2, a1, db)
        return r * float(radius)

    def greatCircle(self, bearing):
        '''Compute the vector normal to great circle obtained by heading
           on the given initial bearing from this point.

           Direction of vector is such that initial bearing vector
           b = c × n, where n is an n-vector representing this point.

           @param bearing: Bearing from this point (compass C{degrees360}).

           @return: Vector representing great circle (L{Vector3d}).

           @example:

           >>> p = LatLon(53.3206, -1.7297)
           >>> g = p.greatCircle(96.0)
           >>> g.toStr()  # (-0.794, 0.129, 0.594)
        '''
        a, b = self.to2ab()
        t = radians(bearing)

        sa, ca, sb, cb, st, ct = sincos2(a, b, t)

        return Vector3d(sb * ct - cb * sa * st,
                       -cb * ct - sb * sa * st,
                        ca * st)  # XXX .unit()?

    def initialBearingTo(self, other, wrap=False, raiser=False):
        '''Compute the initial bearing (forward azimuth) from this
           to an other point.

           @param other: The other point (spherical L{LatLon}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).
           @keyword raiser: Optionally, raise L{CrossError} (C{bool}),
                            use B{C{raiser}}=C{True} for behavior like
                            C{sphericalNvector.LatLon.initialBearingTo}.

           @return: Initial bearing (compass C{degrees360}).

           @raise CrossError: If this and the B{C{other}} point coincide,
                              provided B{C{raiser}} is C{True} and
                              L{crosserrors} is C{True}.

           @raise TypeError: The B{C{other}} point is not L{LatLon}.

           @example:

           >>> p1 = LatLon(52.205, 0.119)
           >>> p2 = LatLon(48.857, 2.351)
           >>> b = p1.initialBearingTo(p2)  # 156.2

           @JSname: I{bearingTo}.
        '''
        self.others(other)

        a1, b1 = self.to2ab()
        a2, b2 = other.to2ab()

        # XXX behavior like sphericalNvector.LatLon.initialBearingTo
        if raiser and crosserrors() and max(abs(a2 - a1), abs(b2 - b1)) < EPS:
            raise CrossError('%s %s: %r' % ('coincident', 'points', other))

        return degrees(bearing_(a1, b1, a2, b2, final=False, wrap=wrap))

    def intermediateTo(self, other, fraction, height=None, wrap=False):
        '''Locate the point at given fraction between this and an
           other point.

           @param other: The other point (L{LatLon}).
           @param fraction: Fraction between both points (float, 0.0 =
                            this point, 1.0 = the other point).
           @keyword height: Optional height, overriding the fractional
                            height (C{meter}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: Intermediate point (L{LatLon}).

           @raise TypeError: The B{C{other}} point is not L{LatLon}.

           @example:

           >>> p1 = LatLon(52.205, 0.119)
           >>> p2 = LatLon(48.857, 2.351)
           >>> p = p1.intermediateTo(p2, 0.25)  # 51.3721°N, 000.7073°E

           @JSname: I{intermediatePointTo}.
        '''
        self.others(other)

        a1, b1 = self.to2ab()
        a2, b2 = other.to2ab()

        db, b2 = unrollPI(b1, b2, wrap=wrap)
        r = haversine_(a2, a1, db)
        sr = sin(r)
        if abs(sr) > EPS:
            sa1, ca1, sa2, ca2, \
            sb1, cb1, sb2, cb2 = sincos2(a1, a2, b1, b2)

            A = sin((1 - fraction) * r) / sr
            B = sin(     fraction  * r) / sr

            x = A * ca1 * cb1 + B * ca2 * cb2
            y = A * ca1 * sb1 + B * ca2 * sb2
            z = A * sa1       + B * sa2

            a = atan2(z, hypot(x, y))
            b = atan2(y, x)

        else:  # points too close
            a = favg(a1, a2, f=fraction)
            b = favg(b1, b2, f=fraction)

        if height is None:
            h = self._havg(other, f=fraction)
        else:
            h = height
        return self.classof(degrees90(a), degrees180(b), height=h)

    def intersection(self, end1, start2, end2,
                           height=None, wrap=False):
        '''Locate the intersection point of two paths both defined
           by two points or a start point and bearing from North.

           @param end1: End point of the first path (L{LatLon}) or
                        the initial bearing at this point (compass
                        C{degrees360}).
           @param start2: Start point of the second path (L{LatLon}).
           @param end2: End point of the second path (L{LatLon}) or
                        the initial bearing at the second point
                        (compass C{degrees360}).
           @keyword height: Optional height for intersection point,
                            overriding the mean height (C{meter}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: The intersection point (L{LatLon}).  An alternate
                    intersection point might be the L{antipode} to
                    the returned result.

           @raise TypeError: The B{C{start2}}, B{C{end1}} or B{C{end2}} is
                             not L{LatLon}.

           @raise ValueError: Intersection is ambiguous or infinite or
                              the paths are parallel, coincident or null.

           @example:

           >>> p = LatLon(51.8853, 0.2545)
           >>> s = LatLon(49.0034, 2.5735)
           >>> i = p.intersection(108.547, s, 32.435)  # '50.9078°N, 004.5084°E'
        '''
        return intersection(self, end1, start2, end2,
                                  height=height, wrap=wrap,
                                  LatLon=self.classof)

    def isenclosedBy(self, points):
        '''Check whether a (convex) polygon encloses this point.

           @param points: The polygon points (L{LatLon}[]).

           @return: C{True} if the polygon encloses this point,
                    C{False} otherwise.

           @raise ValueError: Insufficient number of B{C{points}} or
                              non-convex polygon.

           @raise TypeError: Some B{C{points}} are not L{LatLon}.

           @example:

           >>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
           >>> p = LatLon(45,1, 1.1)
           >>> inside = p.isEnclosedBy(b)  # True
        '''
        n, points = self.points2(points, closed=True)

        n0 = self.toVector3d()

        if iterNumpy2(points):

            v1 = points[n-1].toVector3d()
            v2 = points[n-2].toVector3d()
            gc1 = v2.cross(v1)
            t0 = gc1.angleTo(n0) > PI_2
            for i in range(n):
                v2 = points[i].toVector3d()
                gc = v1.cross(v2)
                v1 = v2

                ti = gc.angleTo(n0) > PI_2
                if ti != t0:
                    return False  # outside

                if gc1.angleTo(gc, vSign=n0) < 0:
                    raise ValueError('non-convex: %r...' % (points[:2],))
                gc1 = gc

        else:
            # get great-circle vector for each edge
            gc, v1 = [], points[n-1].toVector3d()
            for i in range(n):
                v2 = points[i].toVector3d()
                gc.append(v1.cross(v2))
                v1 = v2

            # check whether this point on same side of all
            # polygon edges (to the left or right depending
            # on anti-/clockwise polygon direction)
            t0 = gc[0].angleTo(n0) > PI_2  # True if on the right
            for i in range(1, n):
                ti = gc[i].angleTo(n0) > PI_2
                if ti != t0:  # different sides of edge i
                    return False  # outside

            # check for convex polygon (otherwise
            # the test above is not reliable)
            gc1 = gc[n-1]
            for gc2 in gc:
                # angle between gc vectors, signed by direction of n0
                if gc1.angleTo(gc2, vSign=n0) < 0:
                    raise ValueError('non-convex: %r...' % (points[:2],))
                gc1 = gc2

        return True  # inside

    def isEnclosedBy(self, points):
        '''DEPRECATED, use method C{isenclosedBy}.
        '''
        return self.isenclosedBy(points)

    def midpointTo(self, other, height=None, wrap=False):
        '''Find the midpoint between this and an other point.

           @param other: The other point (L{LatLon}).
           @keyword height: Optional height for midpoint, overriding
                            the mean height (C{meter}).
           @keyword wrap: Wrap and unroll longitudes (C{bool}).

           @return: Midpoint (L{LatLon}).

           @raise TypeError: The B{C{other}} point is not L{LatLon}.

           @example:

           >>> p1 = LatLon(52.205, 0.119)
           >>> p2 = LatLon(48.857, 2.351)
           >>> m = p1.midpointTo(p2)  # '50.5363°N, 001.2746°E'
        '''
        self.others(other)

        # see <https://MathForum.org/library/drmath/view/51822.html>
        a1, b1 = self.to2ab()
        a2, b2 = other.to2ab()

        db, b2 = unrollPI(b1, b2, wrap=wrap)

        sa1, ca1, sa2, ca2, sdb, cdb = sincos2(a1, a2, db)

        x = ca2 * cdb + ca1
        y = ca2 * sdb

        a = atan2(sa1 + sa2, hypot(x, y))
        b = atan2(y, x) + b1

        h = self._havg(other) if height is None else height
        return self.classof(degrees90(a), degrees180(b), height=h)

    def nearestOn(self, point1, point2, radius=R_M, **options):
        '''Locate the point between two points closest and this point.

           Distances are approximated by function L{equirectangular_},
           subject to the supplied B{C{options}}.

           @param point1: Start point (L{LatLon}).
           @param point2: End point (L{LatLon}).
           @keyword radius: Optional, mean earth radius (C{meter}).
           @keyword options: Optional keyword arguments for function
                             L{equirectangular_}.

           @return: Closest point on the arc (L{LatLon}).

           @raise LimitError: Lat- and/or longitudinal delta exceeds
                              B{C{limit}}, see function L{equirectangular_}.

           @raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.

           @see: Functions L{equirectangular_} and L{nearestOn5} and
                 method L{sphericalTrigonometry.LatLon.nearestOn3}.
        '''
        return self.nearestOn3([point1, point2], closed=False, radius=radius,
                                               **options)[0]

    def nearestOn2(self, points, closed=False, radius=R_M, **options):
        '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.

           @return: ... 2-Tuple C{(closest, distance)} of the closest
                    point (L{LatLon}) on the polygon and the distance
                    to that point ...
        '''
        return self.nearestOn3(points, closed=closed, radius=radius,
                                     **options)[:2]

    def nearestOn3(self, points, closed=False, radius=R_M, **options):
        '''Locate the point on a polygon closest to this point.

           Distances are approximated by function L{equirectangular_},
           subject to the supplied B{C{options}}.

           @param points: The polygon points (L{LatLon}[]).
           @keyword closed: Optionally, close the polygon (C{bool}).
           @keyword radius: Optional, mean earth radius (C{meter}).
           @keyword options: Optional keyword arguments for function
                             L{equirectangular_}.

           @return: A L{NearestOn3Tuple}C{(closest, distance, angle)}
                    where C{distance} is the L{equirectangular_} distance
                    between this and the C{closest} point in C{meter},
                    same units as B{C{radius}}.  The C{angle} from this to
                    the C{closest} point is in compass C{degrees360},
                    like function L{compassAngle}.

           @raise LimitError: Lat- and/or longitudinal delta exceeds
                              B{C{limit}}, see function L{equirectangular_}.

           @raise TypeError: Some B{C{points}} are not C{LatLon}.

           @raise ValueError: Insufficient number of B{C{points}}.

           @see: Functions L{compassAngle}, L{equirectangular_} and
                 L{nearestOn5}.
        '''
        a, b, d, c, h = _nearestOn5(self, points, closed=closed, **options)
        return NearestOn3Tuple(self.classof(a, b, height=h),
                               degrees2m(d, radius=radius), c)

    def toCartesian(self, **kwds):  # PYCHOK Cartesian=Cartesian
        '''Convert this point to C{Karney}-based cartesian (ECEF)
           coordinates.

           @keyword kwds: Optional, additional B{C{Cartesian}} keyword
                          arguments, ignored if C{B{Cartesian}=None}.
                          For example, use C{Cartesian=...} to override
                          the L{Cartesian} (sub-)class or specify

           @return: The B{C{Cartesian}} point (L{Cartesian}) or if
                    C{B{Cartesian}=None}, an L{Ecef9Tuple}C{(x, y, z,
                    lat, lon, height, C, M, datum)} with C{C} and C{M}
                    if available.

           @raise TypeError: Invalid B{C{Cartesian}}, B{C{datum}}
                             or B{C{kwds}}.
        '''
        kwds = _2kwds(kwds, Cartesian=Cartesian, datum=self.datum)
        return LatLonSphericalBase.toCartesian(self, **kwds)


_Trll = LatLon(0, 0)  #: (INTERNAL) Reference instance (L{LatLon}).


def _destination2_(a, b, r, t):
    '''(INTERNAL) Computes destination lat- and longitude.

       @param a: Latitude (C{radians}).
       @param b: Longitude (C{radians}).
       @param r: Angular distance (C{radians}).
       @param t: Bearing (compass C{radians}).

       @return: 2-Tuple (lat, lon) of (radians, radians).
    '''
    # see <https://www.EdWilliams.org/avform.htm#LL>
    sa, ca, sr, cr, st, ct = sincos2(a, r, t)

    a = asin(ct * sr * ca + cr * sa)
    d = atan2(st * sr * ca, cr - sa * sin(a))
    # note, in EdWilliams.org/avform.htm W is + and E is -
    return a, b + d


def _destination2(a, b, r, t):
    '''(INTERNAL) Computes destination lat- and longitude.

       @param a: Latitude (C{radians}).
       @param b: Longitude (C{radians}).
       @param r: Angular distance (C{radians}).
       @param t: Bearing (compass C{radians}).

       @return: 2-Tuple (lat, lon) of (C{degrees90}, C{degrees180}).
    '''
    a, b = _destination2_(a, b, r, t)
    return degrees90(a), degrees180(b)


def areaOf(points, radius=R_M, wrap=True):
    '''Calculate the area of a (spherical) polygon (with great circle
       arcs joining the points).

       @param points: The polygon points (L{LatLon}[]).
       @keyword radius: Optional, mean earth radius (C{meter}).
       @keyword wrap: Wrap and unroll longitudes (C{bool}).

       @return: Polygon area (C{meter}, same units as B{C{radius}}, squared).

       @raise TypeError: Some B{C{points}} are not L{LatLon}.

       @raise ValueError: Insufficient number of B{C{points}}.

       @note: The area is based on Karney's U{'Area of a spherical polygon'
              <https://OSGeo-org.1560.x6.nabble.com/
              Area-of-a-spherical-polygon-td3841625.html>}.

       @see: L{pygeodesy.areaOf}, L{sphericalNvector.areaOf} and
             L{ellipsoidalKarney.areaOf}.

       @example:

       >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
       >>> areaOf(b)  # 8666058750.718977

       >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
       >>> areaOf(c)  # 6.18e9
    '''
    n, points = _Trll.points2(points, closed=True)

    # Area method due to Karney: for each edge of the polygon,
    #
    #                tan(Δλ/2) · (tan(φ1/2) + tan(φ2/2))
    #     tan(E/2) = ------------------------------------
    #                     1 + tan(φ1/2) · tan(φ2/2)
    #
    # where E is the spherical excess of the trapezium obtained by
    # extending the edge to the equator-circle vector for each edge

    def _exs(n, points):  # iterate over spherical edge excess
        a1, b1 = points[n-1].to2ab()
        ta1 = tan_2(a1)
        for i in range(n):
            a2, b2 = points[i].to2ab()
            db, b2 = unrollPI(b1, b2, wrap=wrap)
            ta2, tdb = map1(tan_2, a2, db)
            yield atan2(tdb * (ta1 + ta2), 1 + ta1 * ta2)
            ta1, b1 = ta2, b2

    s = fsum(_exs(n, points)) * 2

    if isPoleEnclosedBy(points):
        s = abs(s) - PI2

    return abs(s * radius**2)


def _x3d2(start, end, wrap, n, hs):
    # see <https://www.EdWilliams.org/intersect.htm> (5) ff
    a1, b1 = start.to2ab()

    if isscalar(end):  # bearing, make a point
        a2, b2 = _destination2_(a1, b1, PI_4, radians(end))
    else:  # must be a point
        _Trll.others(end, name='end' + n)
        hs.append(end.height)
        a2, b2 = end.to2ab()

    db, b2 = unrollPI(b1, b2, wrap=wrap)
    if max(abs(db), abs(a2 - a1)) < EPS:
        raise ValueError('intersection %s%s null: %r' % ('path', n, (start, end)))

    # note, in EdWilliams.org/avform.htm W is + and E is -
    b21, b12 = db * 0.5, -(b1 + b2) * 0.5

    sb21, cb21, sb12, cb12, \
    sa21,    _, sa12,    _ = sincos2(b21, b12, a1 - a2, a1 + a2)

    x = Vector3d(sa21 * sb12 * cb21 - sa12 * cb12 * sb21,
                 sa21 * cb12 * cb21 + sa12 * sb12 * sb21,
                 cos(a1) * cos(a2) * sin(db), ll=start)
    return x.unit(), (db, (a2 - a1))  # negated d


def _xb(a1, b1, end, a, b, wrap):
    # difference between the bearing to (a, b) and the given
    # bearing is negative if both are in opposite directions
    r = bearing_(a1, b1, radians(a), radians(b), wrap=wrap)
    return PI_2 - abs(wrapPI(r - radians(end)))


def _xdot(d, a1, b1, a, b, wrap):
    # compute dot product d . (-b + b1, a - a1)
    db, _ = unrollPI(b1, radians(b), wrap=wrap)
    return fdot(d, db, radians(a) - a1)


def intersection(start1, end1, start2, end2,
                 height=None, wrap=False, LatLon=LatLon):
    '''Compute the intersection point of two paths both defined
       by two points or a start point and bearing from North.

       @param start1: Start point of the first path (L{LatLon}).
       @param end1: End point ofthe first path (L{LatLon}) or
                    the initial bearing at the first start point
                    (compass C{degrees360}).
       @param start2: Start point of the second path (L{LatLon}).
       @param end2: End point of the second path (L{LatLon}) or
                    the initial bearing at the second start point
                    (compass C{degrees360}).
       @keyword height: Optional height for the intersection point,
                        overriding the mean height (C{meter}).
       @keyword wrap: Wrap and unroll longitudes (C{bool}).
       @keyword LatLon: Optional (sub-)class to return the intersection
                        point (L{LatLon}) or C{None}.

       @return: The intersection point (B{C{LatLon}}) or a
                L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
                is C{None}.  An alternate intersection point might
                be the L{antipode} to the returned result.

       @raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}.

       @raise ValueError: Intersection is ambiguous or infinite or
                          the paths are parallel, coincident or null.

       @example:

       >>> p = LatLon(51.8853, 0.2545)
       >>> s = LatLon(49.0034, 2.5735)
       >>> i = intersection(p, 108.547, s, 32.435)  # '50.9078°N, 004.5084°E'
    '''
    _Trll.others(start1, name='start1')
    _Trll.others(start2, name='start2')

    hs = [start1.height, start2. height]

    a1, b1 = start1.to2ab()
    a2, b2 = start2.to2ab()

    db, b2 = unrollPI(b1, b2, wrap=wrap)
    r12 = haversine_(a2, a1, db)
    if abs(r12) < EPS:  # [nearly] coincident points
        a, b = map1(degrees, favg(a1, a2), favg(b1, b2))

    # see <https://www.EdWilliams.org/avform.htm#Intersection>
    elif isscalar(end1) and isscalar(end2):  # both bearings
        sa1, ca1, sa2, ca2, sr12, cr12 = sincos2(a1, a2, r12)

        x1, x2 = (sr12 * ca1), (sr12 * ca2)
        if abs(x1) < EPS or abs(x2) < EPS:
            raise ValueError('intersection %s: %r vs %r' % ('parallel',
                             (start1, end1), (start2, end2)))

        # handle domain error for equivalent longitudes,
        # see also functions asin_safe and acos_safe at
        # <https://www.EdWilliams.org/avform.htm#Math>
        t1, t2 = map1(acos1, (sa2 - sa1 * cr12) / x1,
                             (sa1 - sa2 * cr12) / x2)
        if sin(db) > 0:
            t12, t21 = t1, PI2 - t2
        else:
            t12, t21 = PI2 - t1, t2

        t13, t23 = map1(radiansPI2, end1, end2)
        x1, x2 = map1(wrapPI, t13 - t12,  # angle 2-1-3
                              t21 - t23)  # angle 1-2-3
        sx1, cx1, sx2, cx2 = sincos2(x1, x2)
        if sx1 == 0 and sx2 == 0:  # max(abs(sx1), abs(sx2)) < EPS
            raise ValueError('intersection %s: %r vs %r' % ('infinite',
                             (start1, end1), (start2, end2)))
        sx3 = sx1 * sx2
#       if sx3 < 0:
#           raise ValueError('intersection %s: %r vs %r' % ('ambiguous',
#                            (start1, end1), (start2, end2)))
        x3 = acos1(cr12 * sx3 - cx2 * cx1)
        r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))

        a, b = _destination2(a1, b1, r13, t13)
        # choose antipode for opposing bearings
        if _xb(a1, b1, end1, a, b, wrap) < 0 or \
           _xb(a2, b2, end2, a, b, wrap) < 0:
            a, b = antipode(a, b)

    else:  # end point(s) or bearing(s)
        x1, d1 = _x3d2(start1, end1, wrap, '1', hs)
        x2, d2 = _x3d2(start2, end2, wrap, '2', hs)
        x = x1.cross(x2)
        if x.length < EPS:  # [nearly] colinear or parallel paths
            raise ValueError('intersection %s: %r vs %r' % ('colinear',
                             (start1, end1), (start2, end2)))
        a, b = x.to2ll()
        # choose intersection similar to sphericalNvector
        d1 = _xdot(d1, a1, b1, a, b, wrap)
        d2 = _xdot(d2, a2, b2, a, b, wrap)
        if (d1 < 0 and d2 > 0) or (d1 > 0 and d2 < 0):
            a, b = antipode(a, b)

    h = fmean(hs) if height is None else height
    r = LatLon3Tuple(a, b, h) if LatLon is None else \
              LatLon(a, b, height=h)
    return _xnamed(r, intersection.__name__)


def isPoleEnclosedBy(points, wrap=False):
    '''DEPRECATED, use function L{ispolar}.
    '''
    return ispolar(points, wrap=wrap)


def meanOf(points, height=None, LatLon=LatLon):
    '''Compute the geographic mean of several points.

       @param points: Points to be averaged (L{LatLon}[]).
       @keyword height: Optional height at mean point, overriding
                        the mean height (C{meter}).
       @keyword LatLon: Optional (sub-)class to return the mean
                        point (L{LatLon}) or C{None}.

       @return: Point at geographic mean and height (B{C{LatLon}}) or
                a L{LatLon3Tuple}C{(lat, lon, height)} if
                B{C{LatLon}} is C{None}.

       @raise TypeError: Some B{C{points}} are not L{LatLon}.

       @raise ValueError: No B{C{points}}.
    '''
    # geographic mean
    n, points = _Trll.points2(points, closed=False)

    m = sumOf(points[i].toVector3d() for i in range(n))
    a, b = m.to2ll()

    if height is None:
        h = fmean(points[i].height for i in range(n))
    else:
        h = height
    r = LatLon3Tuple(a, b, h) if LatLon is None else \
              LatLon(a, b, height=h)
    return _xnamed(r, meanOf.__name__)


def nearestOn2(point, points, closed=False, radius=R_M,
                              LatLon=LatLon, **options):
    '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.

       @return: ... C{closest} as B{C{LatLon}} or a 2-tuple C{(lat, lon)}
                without the height if B{C{LatLon}} is C{None} ...
    '''
    a, b, d, _, h = _nearestOn5(point, points, closed=closed, **options)
    ll = (a, b) if LatLon is None else LatLon(a, b, height=h)
    return ll, degrees2m(d, radius=radius)


def nearestOn3(point, points, closed=False, radius=R_M,
                              LatLon=LatLon, **options):
    '''Locate the point on a polygon closest to an other, reference point.

       Distances are approximated by function L{equirectangular_},
       subject to the supplied B{C{options}}.

       @param point: The other, reference point (L{LatLon}).
       @param points: The polygon points (L{LatLon}[]).
       @keyword closed: Optionally, close the polygon (C{bool}).
       @keyword radius: Optional, mean earth radius (C{meter}).
       @keyword LatLon: Optional (sub-)class to return the closest
                        point (L{LatLon}) or C{None}.
       @keyword options: Optional keyword arguments for function
                         L{equirectangular_}.

       @return: A L{NearestOn3Tuple}C{(closest, distance, angle)}.  The
                C{distance} is the L{equirectangular_} distance between
                the C{closest} and reference B{C{point}} in C{meter}, same
                units as B{C{radius}}.  The C{angle} from the reference
                B{C{point}} to the C{closest} is in compass C{degrees360},
                like function L{compassAngle}.  The C{height} is the
                (interpolated) height at the C{closest} point.

       @raise LimitError: Lat- and/or longitudinal delta exceeds the
                          B{C{limit}}, see function L{equirectangular_}.

       @raise TypeError: Some I{points} are not C{LatLon}.

       @raise ValueError: Insufficient number of B{C{points}}.

       @see: Functions L{equirectangular_} and L{nearestOn5}.
    '''
    a, b, d, c, h = _nearestOn5(point, points, closed=closed,
                                               LatLon=None, **options)
    r = LatLon3Tuple(a, b, h) if LatLon is None else \
              LatLon(a, b, height=h)
    r = NearestOn3Tuple(r, degrees2m(d, radius=radius), c)
    return _xnamed(r, nearestOn3.__name__)


def perimeterOf(points, closed=False, radius=R_M, wrap=True):
    '''Compute the perimeter of a (spherical) polygon (with great circle
       arcs joining the points).

       @param points: The polygon points (L{LatLon}[]).
       @keyword closed: Optionally, close the polygon (C{bool}).
       @keyword radius: Optional, mean earth radius (C{meter}).
       @keyword wrap: Wrap and unroll longitudes (C{bool}).

       @return: Polygon perimeter (C{meter}, same units as B{C{radius}}).

       @raise TypeError: Some B{C{points}} are not L{LatLon}.

       @raise ValueError: Insufficient number of B{C{points}}.

       @note: This perimeter is based on the L{haversine} formula.

       @see: L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
             and L{ellipsoidalKarney.perimeterOf}.
    '''
    n, points = _Trll.points2(points, closed=closed)

    def _rads(n, points, closed):  # angular edge lengths in radians
        i, m = _imdex2(closed, n)
        a1, b1 = points[i].to2ab()
        for i in range(m, n):
            a2, b2 = points[i].to2ab()
            db, b2 = unrollPI(b1, b2, wrap=wrap)
            yield haversine_(a2, a1, db)
            a1, b1 = a2, b2

    r = fsum(_rads(n, points, closed))
    return r * float(radius)

# **) MIT License
#
# Copyright (C) 2016-2020 -- mrJean1 at Gmail -- All Rights Reserved.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.
