fix subject f string

streamlit-venv/lib/python3.10/site-packages/geographiclib/__init__.py

@@ -0,0 +1,12 @@
+"""geographiclib: geodesic routines from GeographicLib"""
+
+__version_info__ = (2,
+                    0,
+                    0)
+"""GeographicLib version as a tuple"""
+
+__version__ = "2.0"
+"""GeographicLib version as a string"""
+
+__release_date__ = "2022-04-23"
+"""GeographicLib release date"""

streamlit-venv/lib/python3.10/site-packages/geographiclib/accumulator.py

@@ -0,0 +1,87 @@
+"""accumulator.py: transcription of GeographicLib::Accumulator class."""
+# accumulator.py
+#
+# This is a rather literal translation of the GeographicLib::Accumulator class
+# from to python.  See the documentation for the C++ class for more information
+# at
+#
+#    https://geographiclib.sourceforge.io/html/annotated.html
+#
+# Copyright (c) Charles Karney (2011-2019) <charles@karney.com> and
+# licensed under the MIT/X11 License.  For more information, see
+# https://geographiclib.sourceforge.io/
+######################################################################
+
+from geographiclib.geomath import Math
+
+class Accumulator:
+  """Like math.fsum, but allows a running sum"""
+
+  def Set(self, y):
+    """Set value from argument"""
+    if isinstance(y, Accumulator):
+      self._s, self._t = y._s, y._t
+    else:
+      self._s, self._t = float(y), 0.0
+
+  def __init__(self, y = 0.0):
+    """Constructor"""
+    self._s = self._t = 0.0
+    self.Set(y)
+
+  def Add(self, y):
+    """Add a value"""
+    # Here's Shewchuk's solution...
+    # hold exact sum as [s, t, u]
+    y, u = Math.sum(y, self._t)             # Accumulate starting at
+    self._s, self._t = Math.sum(y, self._s) # least significant end
+    # Start is _s, _t decreasing and non-adjacent.  Sum is now (s + t + u)
+    # exactly with s, t, u non-adjacent and in decreasing order (except
+    # for possible zeros).  The following code tries to normalize the
+    # result.  Ideally, we want _s = round(s+t+u) and _u = round(s+t+u -
+    # _s).  The follow does an approximate job (and maintains the
+    # decreasing non-adjacent property).  Here are two "failures" using
+    # 3-bit floats:
+    #
+    # Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
+    # [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
+    #
+    # Case 2: _s+_t is not as close to s+t+u as it shold be
+    # [64, 5] + 4 -> [64, 8, 1] -> [64,  8] = 72 (off by 1)
+    #                    should be [80, -7] = 73 (exact)
+    #
+    # "Fixing" these problems is probably not worth the expense.  The
+    # representation inevitably leads to small errors in the accumulated
+    # values.  The additional errors illustrated here amount to 1 ulp of
+    # the less significant word during each addition to the Accumulator
+    # and an additional possible error of 1 ulp in the reported sum.
+    #
+    # Incidentally, the "ideal" representation described above is not
+    # canonical, because _s = round(_s + _t) may not be true.  For
+    # example, with 3-bit floats:
+    #
+    # [128, 16] + 1 -> [160, -16] -- 160 = round(145).
+    # But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
+    #
+    if self._s == 0:            # This implies t == 0,
+      self._s = u               # so result is u
+    else:
+      self._t += u              # otherwise just accumulate u to t.
+
+  def Sum(self, y = 0.0):
+    """Return sum + y"""
+    if y == 0.0:
+      return self._s
+    b = Accumulator(self)
+    b.Add(y)
+    return b._s
+
+  def Negate(self):
+    """Negate sum"""
+    self._s *= -1
+    self._t *= -1
+
+  def Remainder(self, y):
+    """Remainder on division by y"""
+    self._s = Math.remainder(self._s, y)
+    self.Add(0.0)

streamlit-venv/lib/python3.10/site-packages/geographiclib/constants.py

@@ -0,0 +1,22 @@
+"""Define the WGS84 ellipsoid"""
+# constants.py
+#
+# This is a translation of the GeographicLib::Constants class to python.  See
+# the documentation for the C++ class for more information at
+#
+#    https://geographiclib.sourceforge.io/html/annotated.html
+#
+# Copyright (c) Charles Karney (2011-2016) <charles@karney.com> and
+# licensed under the MIT/X11 License.  For more information, see
+# https://geographiclib.sourceforge.io/
+######################################################################
+
+class Constants:
+  """
+  Constants describing the WGS84 ellipsoid
+  """
+
+  WGS84_a = 6378137.0           # meters
+  """the equatorial radius in meters of the WGS84 ellipsoid in meters"""
+  WGS84_f = 1/298.257223563
+  """the flattening of the WGS84 ellipsoid, 1/298.257223563"""

streamlit-venv/lib/python3.10/site-packages/geographiclib/geodesiccapability.py

@@ -0,0 +1,41 @@
+"""geodesiccapability.py: capability constants for geodesic{,line}.py"""
+# geodesiccapability.py
+#
+# This gathers the capability constants need by geodesic.py and
+# geodesicline.py.  See the documentation for the GeographicLib::Geodesic class
+# for more information at
+#
+#    https://geographiclib.sourceforge.io/html/annotated.html
+#
+# Copyright (c) Charles Karney (2011-2014) <charles@karney.com> and licensed
+# under the MIT/X11 License.  For more information, see
+# https://geographiclib.sourceforge.io/
+######################################################################
+
+class GeodesicCapability:
+  """
+  Capability constants shared between Geodesic and GeodesicLine.
+  """
+
+  CAP_NONE = 0
+  CAP_C1   = 1 << 0
+  CAP_C1p  = 1 << 1
+  CAP_C2   = 1 << 2
+  CAP_C3   = 1 << 3
+  CAP_C4   = 1 << 4
+  CAP_ALL  = 0x1F
+  CAP_MASK = CAP_ALL
+  OUT_ALL  = 0x7F80
+  OUT_MASK = 0xFF80             # Includes LONG_UNROLL
+  EMPTY         = 0
+  LATITUDE      = 1 << 7  | CAP_NONE
+  LONGITUDE     = 1 << 8  | CAP_C3
+  AZIMUTH       = 1 << 9  | CAP_NONE
+  DISTANCE      = 1 << 10 | CAP_C1
+  STANDARD      = LATITUDE | LONGITUDE | AZIMUTH | DISTANCE
+  DISTANCE_IN   = 1 << 11 | CAP_C1 | CAP_C1p
+  REDUCEDLENGTH = 1 << 12 | CAP_C1 | CAP_C2
+  GEODESICSCALE = 1 << 13 | CAP_C1 | CAP_C2
+  AREA          = 1 << 14 | CAP_C4
+  LONG_UNROLL   = 1 << 15
+  ALL           = OUT_ALL | CAP_ALL # Does not include LONG_UNROLL

streamlit-venv/lib/python3.10/site-packages/geographiclib/geodesicline.py

@@ -0,0 +1,427 @@
+"""Define the :class:`~geographiclib.geodesicline.GeodesicLine` class
+
+The constructor defines the starting point of the line.  Points on the
+line are given by
+
+  * :meth:`~geographiclib.geodesicline.GeodesicLine.Position` position
+    given in terms of distance
+  * :meth:`~geographiclib.geodesicline.GeodesicLine.ArcPosition` position
+    given in terms of spherical arc length
+
+A reference point 3 can be defined with
+
+  * :meth:`~geographiclib.geodesicline.GeodesicLine.SetDistance` set
+    position of 3 in terms of the distance from the starting point
+  * :meth:`~geographiclib.geodesicline.GeodesicLine.SetArc` set
+    position of 3 in terms of the spherical arc length from the starting point
+
+The object can also be constructed by
+
+  * :meth:`Geodesic.Line <geographiclib.geodesic.Geodesic.Line>`
+  * :meth:`Geodesic.DirectLine <geographiclib.geodesic.Geodesic.DirectLine>`
+  * :meth:`Geodesic.ArcDirectLine
+    <geographiclib.geodesic.Geodesic.ArcDirectLine>`
+  * :meth:`Geodesic.InverseLine <geographiclib.geodesic.Geodesic.InverseLine>`
+
+The public attributes for this class are
+
+  * :attr:`~geographiclib.geodesicline.GeodesicLine.a`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.f`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.caps`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.lat1`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.lon1`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.azi1`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.salp1`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.calp1`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.s13`
+    :attr:`~geographiclib.geodesicline.GeodesicLine.a13`
+
+"""
+# geodesicline.py
+#
+# This is a rather literal translation of the GeographicLib::GeodesicLine class
+# to python.  See the documentation for the C++ class for more information at
+#
+#    https://geographiclib.sourceforge.io/html/annotated.html
+#
+# The algorithms are derived in
+#
+#    Charles F. F. Karney,
+#    Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
+#    https://doi.org/10.1007/s00190-012-0578-z
+#    Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
+#
+# Copyright (c) Charles Karney (2011-2022) <charles@karney.com> and licensed
+# under the MIT/X11 License.  For more information, see
+# https://geographiclib.sourceforge.io/
+######################################################################
+
+import math
+from geographiclib.geomath import Math
+from geographiclib.geodesiccapability import GeodesicCapability
+
+class GeodesicLine:
+  """Points on a geodesic path"""
+
+  def __init__(self, geod, lat1, lon1, azi1,
+               caps = GeodesicCapability.STANDARD |
+               GeodesicCapability.DISTANCE_IN,
+               salp1 = math.nan, calp1 = math.nan):
+    """Construct a GeodesicLine object
+
+    :param geod: a :class:`~geographiclib.geodesic.Geodesic` object
+    :param lat1: latitude of the first point in degrees
+    :param lon1: longitude of the first point in degrees
+    :param azi1: azimuth at the first point in degrees
+    :param caps: the :ref:`capabilities <outmask>`
+
+    This creates an object allowing points along a geodesic starting at
+    (*lat1*, *lon1*), with azimuth *azi1* to be found.  The default
+    value of *caps* is STANDARD | DISTANCE_IN.  The optional parameters
+    *salp1* and *calp1* should not be supplied; they are part of the
+    private interface.
+
+    """
+
+    from geographiclib.geodesic import Geodesic
+    self.a = geod.a
+    """The equatorial radius in meters (readonly)"""
+    self.f = geod.f
+    """The flattening (readonly)"""
+    self._b = geod._b
+    self._c2 = geod._c2
+    self._f1 = geod._f1
+    self.caps = (caps | Geodesic.LATITUDE | Geodesic.AZIMUTH |
+                  Geodesic.LONG_UNROLL)
+    """the capabilities (readonly)"""
+
+    # Guard against underflow in salp0
+    self.lat1 = Math.LatFix(lat1)
+    """the latitude of the first point in degrees (readonly)"""
+    self.lon1 = lon1
+    """the longitude of the first point in degrees (readonly)"""
+    if math.isnan(salp1) or math.isnan(calp1):
+      self.azi1 = Math.AngNormalize(azi1)
+      self.salp1, self.calp1 = Math.sincosd(Math.AngRound(azi1))
+    else:
+      self.azi1 = azi1
+      """the azimuth at the first point in degrees (readonly)"""
+      self.salp1 = salp1
+      """the sine of the azimuth at the first point (readonly)"""
+      self.calp1 = calp1
+      """the cosine of the azimuth at the first point (readonly)"""
+
+    # real cbet1, sbet1
+    sbet1, cbet1 = Math.sincosd(Math.AngRound(self.lat1)); sbet1 *= self._f1
+    # Ensure cbet1 = +epsilon at poles
+    sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
+    self._dn1 = math.sqrt(1 + geod._ep2 * Math.sq(sbet1))
+
+    # Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
+    self._salp0 = self.salp1 * cbet1 # alp0 in [0, pi/2 - |bet1|]
+    # Alt: calp0 = hypot(sbet1, calp1 * cbet1).  The following
+    # is slightly better (consider the case salp1 = 0).
+    self._calp0 = math.hypot(self.calp1, self.salp1 * sbet1)
+    # Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
+    # sig = 0 is nearest northward crossing of equator.
+    # With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
+    # With bet1 =  pi/2, alp1 = -pi, sig1 =  pi/2
+    # With bet1 = -pi/2, alp1 =  0 , sig1 = -pi/2
+    # Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
+    # With alp0 in (0, pi/2], quadrants for sig and omg coincide.
+    # No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
+    # With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
+    self._ssig1 = sbet1; self._somg1 = self._salp0 * sbet1
+    self._csig1 = self._comg1 = (cbet1 * self.calp1
+                                 if sbet1 != 0 or self.calp1 != 0 else 1)
+    # sig1 in (-pi, pi]
+    self._ssig1, self._csig1 = Math.norm(self._ssig1, self._csig1)
+    # No need to normalize
+    # self._somg1, self._comg1 = Math.norm(self._somg1, self._comg1)
+
+    self._k2 = Math.sq(self._calp0) * geod._ep2
+    eps = self._k2 / (2 * (1 + math.sqrt(1 + self._k2)) + self._k2)
+
+    if self.caps & Geodesic.CAP_C1:
+      self._A1m1 = Geodesic._A1m1f(eps)
+      self._C1a = list(range(Geodesic.nC1_ + 1))
+      Geodesic._C1f(eps, self._C1a)
+      self._B11 = Geodesic._SinCosSeries(
+        True, self._ssig1, self._csig1, self._C1a)
+      s = math.sin(self._B11); c = math.cos(self._B11)
+      # tau1 = sig1 + B11
+      self._stau1 = self._ssig1 * c + self._csig1 * s
+      self._ctau1 = self._csig1 * c - self._ssig1 * s
+      # Not necessary because C1pa reverts C1a
+      #    _B11 = -_SinCosSeries(true, _stau1, _ctau1, _C1pa)
+
+    if self.caps & Geodesic.CAP_C1p:
+      self._C1pa = list(range(Geodesic.nC1p_ + 1))
+      Geodesic._C1pf(eps, self._C1pa)
+
+    if self.caps & Geodesic.CAP_C2:
+      self._A2m1 = Geodesic._A2m1f(eps)
+      self._C2a = list(range(Geodesic.nC2_ + 1))
+      Geodesic._C2f(eps, self._C2a)
+      self._B21 = Geodesic._SinCosSeries(
+        True, self._ssig1, self._csig1, self._C2a)
+
+    if self.caps & Geodesic.CAP_C3:
+      self._C3a = list(range(Geodesic.nC3_))
+      geod._C3f(eps, self._C3a)
+      self._A3c = -self.f * self._salp0 * geod._A3f(eps)
+      self._B31 = Geodesic._SinCosSeries(
+        True, self._ssig1, self._csig1, self._C3a)
+
+    if self.caps & Geodesic.CAP_C4:
+      self._C4a = list(range(Geodesic.nC4_))
+      geod._C4f(eps, self._C4a)
+      # Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
+      self._A4 = Math.sq(self.a) * self._calp0 * self._salp0 * geod._e2
+      self._B41 = Geodesic._SinCosSeries(
+        False, self._ssig1, self._csig1, self._C4a)
+    self.s13 = math.nan
+    """the distance between point 1 and point 3 in meters (readonly)"""
+    self.a13 = math.nan
+    """the arc length between point 1 and point 3 in degrees (readonly)"""
+
+  # return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
+  def _GenPosition(self, arcmode, s12_a12, outmask):
+    """Private: General solution of position along geodesic"""
+    from geographiclib.geodesic import Geodesic
+    a12 = lat2 = lon2 = azi2 = s12 = m12 = M12 = M21 = S12 = math.nan
+    outmask &= self.caps & Geodesic.OUT_MASK
+    if not (arcmode or
+            (self.caps & (Geodesic.OUT_MASK & Geodesic.DISTANCE_IN))):
+      # Uninitialized or impossible distance calculation requested
+      return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
+
+    # Avoid warning about uninitialized B12.
+    B12 = 0.0; AB1 = 0.0
+    if arcmode:
+      # Interpret s12_a12 as spherical arc length
+      sig12 = math.radians(s12_a12)
+      ssig12, csig12 = Math.sincosd(s12_a12)
+    else:
+      # Interpret s12_a12 as distance
+      tau12 = s12_a12 / (self._b * (1 + self._A1m1))
+      tau12 = tau12 if math.isfinite(tau12) else math.nan
+      s = math.sin(tau12); c = math.cos(tau12)
+      # tau2 = tau1 + tau12
+      B12 = - Geodesic._SinCosSeries(True,
+                                    self._stau1 * c + self._ctau1 * s,
+                                    self._ctau1 * c - self._stau1 * s,
+                                    self._C1pa)
+      sig12 = tau12 - (B12 - self._B11)
+      ssig12 = math.sin(sig12); csig12 = math.cos(sig12)
+      if abs(self.f) > 0.01:
+        # Reverted distance series is inaccurate for |f| > 1/100, so correct
+        # sig12 with 1 Newton iteration.  The following table shows the
+        # approximate maximum error for a = WGS_a() and various f relative to
+        # GeodesicExact.
+        #     erri = the error in the inverse solution (nm)
+        #     errd = the error in the direct solution (series only) (nm)
+        #     errda = the error in the direct solution (series + 1 Newton) (nm)
+        #
+        #       f     erri  errd errda
+        #     -1/5    12e6 1.2e9  69e6
+        #     -1/10  123e3  12e6 765e3
+        #     -1/20   1110 108e3  7155
+        #     -1/50  18.63 200.9 27.12
+        #     -1/100 18.63 23.78 23.37
+        #     -1/150 18.63 21.05 20.26
+        #      1/150 22.35 24.73 25.83
+        #      1/100 22.35 25.03 25.31
+        #      1/50  29.80 231.9 30.44
+        #      1/20   5376 146e3  10e3
+        #      1/10  829e3  22e6 1.5e6
+        #      1/5   157e6 3.8e9 280e6
+        ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12
+        csig2 = self._csig1 * csig12 - self._ssig1 * ssig12
+        B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
+        serr = ((1 + self._A1m1) * (sig12 + (B12 - self._B11)) -
+                s12_a12 / self._b)
+        sig12 = sig12 - serr / math.sqrt(1 + self._k2 * Math.sq(ssig2))
+        ssig12 = math.sin(sig12); csig12 = math.cos(sig12)
+        # Update B12 below
+
+    # real omg12, lam12, lon12
+    # real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2
+    # sig2 = sig1 + sig12
+    ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12
+    csig2 = self._csig1 * csig12 - self._ssig1 * ssig12
+    dn2 = math.sqrt(1 + self._k2 * Math.sq(ssig2))
+    if outmask & (
+      Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
+      if arcmode or abs(self.f) > 0.01:
+        B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
+      AB1 = (1 + self._A1m1) * (B12 - self._B11)
+    # sin(bet2) = cos(alp0) * sin(sig2)
+    sbet2 = self._calp0 * ssig2
+    # Alt: cbet2 = hypot(csig2, salp0 * ssig2)
+    cbet2 = math.hypot(self._salp0, self._calp0 * csig2)
+    if cbet2 == 0:
+      # I.e., salp0 = 0, csig2 = 0.  Break the degeneracy in this case
+      cbet2 = csig2 = Geodesic.tiny_
+    # tan(alp0) = cos(sig2)*tan(alp2)
+    salp2 = self._salp0; calp2 = self._calp0 * csig2 # No need to normalize
+
+    if outmask & Geodesic.DISTANCE:
+      s12 = self._b * ((1 + self._A1m1) * sig12 + AB1) if arcmode else s12_a12
+
+    if outmask & Geodesic.LONGITUDE:
+      # tan(omg2) = sin(alp0) * tan(sig2)
+      somg2 = self._salp0 * ssig2; comg2 = csig2 # No need to normalize
+      E = math.copysign(1, self._salp0)          # East or west going?
+      # omg12 = omg2 - omg1
+      omg12 = (E * (sig12
+                    - (math.atan2(          ssig2,       csig2) -
+                       math.atan2(    self._ssig1, self._csig1))
+                    + (math.atan2(E *       somg2,       comg2) -
+                       math.atan2(E * self._somg1, self._comg1)))
+               if outmask & Geodesic.LONG_UNROLL
+               else math.atan2(somg2 * self._comg1 - comg2 * self._somg1,
+                               comg2 * self._comg1 + somg2 * self._somg1))
+      lam12 = omg12 + self._A3c * (
+        sig12 + (Geodesic._SinCosSeries(True, ssig2, csig2, self._C3a)
+                 - self._B31))
+      lon12 = math.degrees(lam12)
+      lon2 = (self.lon1 + lon12 if outmask & Geodesic.LONG_UNROLL else
+              Math.AngNormalize(Math.AngNormalize(self.lon1) +
+                                Math.AngNormalize(lon12)))
+
+    if outmask & Geodesic.LATITUDE:
+      lat2 = Math.atan2d(sbet2, self._f1 * cbet2)
+
+    if outmask & Geodesic.AZIMUTH:
+      azi2 = Math.atan2d(salp2, calp2)
+
+    if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
+      B22 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C2a)
+      AB2 = (1 + self._A2m1) * (B22 - self._B21)
+      J12 = (self._A1m1 - self._A2m1) * sig12 + (AB1 - AB2)
+      if outmask & Geodesic.REDUCEDLENGTH:
+        # Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
+        # accurate cancellation in the case of coincident points.
+        m12 = self._b * ((      dn2 * (self._csig1 * ssig2) -
+                          self._dn1 * (self._ssig1 * csig2))
+                         - self._csig1 * csig2 * J12)
+      if outmask & Geodesic.GEODESICSCALE:
+        t = (self._k2 * (ssig2 - self._ssig1) *
+             (ssig2 + self._ssig1) / (self._dn1 + dn2))
+        M12 = csig12 + (t * ssig2 - csig2 * J12) * self._ssig1 / self._dn1
+        M21 = csig12 - (t * self._ssig1 - self._csig1 * J12) * ssig2 / dn2
+
+    if outmask & Geodesic.AREA:
+      B42 = Geodesic._SinCosSeries(False, ssig2, csig2, self._C4a)
+      # real salp12, calp12
+      if self._calp0 == 0 or self._salp0 == 0:
+        # alp12 = alp2 - alp1, used in atan2 so no need to normalize
+        salp12 = salp2 * self.calp1 - calp2 * self.salp1
+        calp12 = calp2 * self.calp1 + salp2 * self.salp1
+      else:
+        # tan(alp) = tan(alp0) * sec(sig)
+        # tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
+        # = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
+        # If csig12 > 0, write
+        #   csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
+        # else
+        #   csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
+        # No need to normalize
+        salp12 = self._calp0 * self._salp0 * (
+          self._csig1 * (1 - csig12) + ssig12 * self._ssig1 if csig12 <= 0
+          else ssig12 * (self._csig1 * ssig12 / (1 + csig12) + self._ssig1))
+        calp12 = (Math.sq(self._salp0) +
+                  Math.sq(self._calp0) * self._csig1 * csig2)
+      S12 = (self._c2 * math.atan2(salp12, calp12) +
+             self._A4 * (B42 - self._B41))
+
+    a12 = s12_a12 if arcmode else math.degrees(sig12)
+    return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
+
+  def Position(self, s12, outmask = GeodesicCapability.STANDARD):
+    """Find the position on the line given *s12*
+
+    :param s12: the distance from the first point to the second in
+      meters
+    :param outmask: the :ref:`output mask <outmask>`
+    :return: a :ref:`dict`
+
+    The default value of *outmask* is STANDARD, i.e., the *lat1*,
+    *lon1*, *azi1*, *lat2*, *lon2*, *azi2*, *s12*, *a12* entries are
+    returned.  The :class:`~geographiclib.geodesicline.GeodesicLine`
+    object must have been constructed with the DISTANCE_IN capability.
+
+    """
+
+    from geographiclib.geodesic import Geodesic
+    result = {'lat1': self.lat1,
+              'lon1': self.lon1 if outmask & Geodesic.LONG_UNROLL else
+              Math.AngNormalize(self.lon1),
+              'azi1': self.azi1, 's12': s12}
+    a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = self._GenPosition(
+      False, s12, outmask)
+    outmask &= Geodesic.OUT_MASK
+    result['a12'] = a12
+    if outmask & Geodesic.LATITUDE: result['lat2'] = lat2
+    if outmask & Geodesic.LONGITUDE: result['lon2'] = lon2
+    if outmask & Geodesic.AZIMUTH: result['azi2'] = azi2
+    if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
+    if outmask & Geodesic.GEODESICSCALE:
+      result['M12'] = M12; result['M21'] = M21
+    if outmask & Geodesic.AREA: result['S12'] = S12
+    return result
+
+  def ArcPosition(self, a12, outmask = GeodesicCapability.STANDARD):
+    """Find the position on the line given *a12*
+
+    :param a12: spherical arc length from the first point to the second
+      in degrees
+    :param outmask: the :ref:`output mask <outmask>`
+    :return: a :ref:`dict`
+
+    The default value of *outmask* is STANDARD, i.e., the *lat1*,
+    *lon1*, *azi1*, *lat2*, *lon2*, *azi2*, *s12*, *a12* entries are
+    returned.
+
+    """
+
+    from geographiclib.geodesic import Geodesic
+    result = {'lat1': self.lat1,
+              'lon1': self.lon1 if outmask & Geodesic.LONG_UNROLL else
+              Math.AngNormalize(self.lon1),
+              'azi1': self.azi1, 'a12': a12}
+    a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = self._GenPosition(
+      True, a12, outmask)
+    outmask &= Geodesic.OUT_MASK
+    if outmask & Geodesic.DISTANCE: result['s12'] = s12
+    if outmask & Geodesic.LATITUDE: result['lat2'] = lat2
+    if outmask & Geodesic.LONGITUDE: result['lon2'] = lon2
+    if outmask & Geodesic.AZIMUTH: result['azi2'] = azi2
+    if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
+    if outmask & Geodesic.GEODESICSCALE:
+      result['M12'] = M12; result['M21'] = M21
+    if outmask & Geodesic.AREA: result['S12'] = S12
+    return result
+
+  def SetDistance(self, s13):
+    """Specify the position of point 3 in terms of distance
+
+    :param s13: distance from point 1 to point 3 in meters
+
+    """
+
+    self.s13 = s13
+    self.a13, _, _, _, _, _, _, _, _ = self._GenPosition(False, self.s13, 0)
+
+  def SetArc(self, a13):
+    """Specify the position of point 3 in terms of arc length
+
+    :param a13: spherical arc length from point 1 to point 3 in degrees
+
+    """
+
+    from geographiclib.geodesic import Geodesic
+    self.a13 = a13
+    _, _, _, _, self.s13, _, _, _, _ = self._GenPosition(True, self.a13,
+                                                         Geodesic.DISTANCE)

streamlit-venv/lib/python3.10/site-packages/geographiclib/geomath.py

@@ -0,0 +1,162 @@
+"""geomath.py: transcription of GeographicLib::Math class."""
+# geomath.py
+#
+# This is a rather literal translation of the GeographicLib::Math class to
+# python.  See the documentation for the C++ class for more information at
+#
+#    https://geographiclib.sourceforge.io/html/annotated.html
+#
+# Copyright (c) Charles Karney (2011-2021) <charles@karney.com> and
+# licensed under the MIT/X11 License.  For more information, see
+# https://geographiclib.sourceforge.io/
+######################################################################
+
+import sys
+import math
+
+class Math:
+  """
+  Additional math routines for GeographicLib.
+  """
+
+  @staticmethod
+  def sq(x):
+    """Square a number"""
+
+    return x * x
+
+  @staticmethod
+  def cbrt(x):
+    """Real cube root of a number"""
+
+    return math.copysign(math.pow(abs(x), 1/3.0), x)
+
+  @staticmethod
+  def norm(x, y):
+    """Private: Normalize a two-vector."""
+
+    r = (math.sqrt(Math.sq(x) + Math.sq(y))
+         # hypot is inaccurate for 3.[89].  Problem reported by agdhruv
+         # https://github.com/geopy/geopy/issues/466 ; see
+         # https://bugs.python.org/issue43088
+         # Visual Studio 2015 32-bit has a similar problem.
+         if (3, 8) <= sys.version_info < (3, 10)
+         else math.hypot(x, y))
+    return x/r, y/r
+
+  @staticmethod
+  def sum(u, v):
+    """Error free transformation of a sum."""
+
+    # Error free transformation of a sum.  Note that t can be the same as one
+    # of the first two arguments.
+    s = u + v
+    up = s - v
+    vpp = s - up
+    up -= u
+    vpp -= v
+    t = s if s == 0 else 0.0 - (up + vpp)
+    # u + v =       s      + t
+    #       = round(u + v) + t
+    return s, t
+
+  @staticmethod
+  def polyval(N, p, s, x):
+    """Evaluate a polynomial."""
+
+    y = float(0 if N < 0 else p[s]) # make sure the returned value is a float
+    while N > 0:
+      N -= 1; s += 1
+      y = y * x + p[s]
+    return y
+
+  @staticmethod
+  def AngRound(x):
+    """Private: Round an angle so that small values underflow to zero."""
+
+    # The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
+    # for reals = 0.7 pm on the earth if x is an angle in degrees.  (This
+    # is about 1000 times more resolution than we get with angles around 90
+    # degrees.)  We use this to avoid having to deal with near singular
+    # cases when x is non-zero but tiny (e.g., 1.0e-200).
+    z = 1/16.0
+    y = abs(x)
+    # The compiler mustn't "simplify" z - (z - y) to y
+    if y < z: y = z - (z - y)
+    return math.copysign(y, x)
+
+  @staticmethod
+  def remainder(x, y):
+    """remainder of x/y in the range [-y/2, y/2]."""
+
+    return math.remainder(x, y) if math.isfinite(x) else math.nan
+
+  @staticmethod
+  def AngNormalize(x):
+    """reduce angle to [-180,180]"""
+
+    y = Math.remainder(x, 360)
+    return math.copysign(180.0, x) if abs(y) == 180 else y
+
+  @staticmethod
+  def LatFix(x):
+    """replace angles outside [-90,90] by NaN"""
+
+    return math.nan if abs(x) > 90 else x
+
+  @staticmethod
+  def AngDiff(x, y):
+    """compute y - x and reduce to [-180,180] accurately"""
+
+    d, t = Math.sum(Math.remainder(-x, 360), Math.remainder(y, 360))
+    d, t = Math.sum(Math.remainder(d, 360), t)
+    if d == 0 or abs(d) == 180:
+      d = math.copysign(d, y - x if t == 0 else -t)
+    return d, t
+
+  @staticmethod
+  def sincosd(x):
+    """Compute sine and cosine of x in degrees."""
+
+    r = math.fmod(x, 360) if math.isfinite(x) else math.nan
+    q = 0 if math.isnan(r) else int(round(r / 90))
+    r -= 90 * q; r = math.radians(r)
+    s = math.sin(r); c = math.cos(r)
+    q = q % 4
+    if   q == 1: s, c =  c, -s
+    elif q == 2: s, c = -s, -c
+    elif q == 3: s, c = -c,  s
+    c = c + 0.0
+    if s == 0: s = math.copysign(s, x)
+    return s, c
+
+  @staticmethod
+  def sincosde(x, t):
+    """Compute sine and cosine of (x + t) in degrees with x in [-180, 180]"""
+
+    q = int(round(x / 90)) if math.isfinite(x) else 0
+    r = x - 90 * q; r = math.radians(Math.AngRound(r + t))
+    s = math.sin(r); c = math.cos(r)
+    q = q % 4
+    if   q == 1: s, c =  c, -s
+    elif q == 2: s, c = -s, -c
+    elif q == 3: s, c = -c,  s
+    c = c + 0.0
+    if s == 0: s = math.copysign(s, x)
+    return s, c
+
+  @staticmethod
+  def atan2d(y, x):
+    """compute atan2(y, x) with the result in degrees"""
+
+    if abs(y) > abs(x):
+      q = 2; x, y = y, x
+    else:
+      q = 0
+    if x < 0:
+      q += 1; x = -x
+    ang = math.degrees(math.atan2(y, x))
+    if   q == 1: ang = math.copysign(180, y) - ang
+    elif q == 2: ang =                90      - ang
+    elif q == 3: ang =               -90     + ang
+    return ang

streamlit-venv/lib/python3.10/site-packages/geographiclib/polygonarea.py

@@ -0,0 +1,326 @@
+"""Define the :class:`~geographiclib.polygonarea.PolygonArea` class
+
+The constructor initializes a empty polygon.  The available methods are
+
+  * :meth:`~geographiclib.polygonarea.PolygonArea.Clear` reset the
+    polygon
+  * :meth:`~geographiclib.polygonarea.PolygonArea.AddPoint` add a vertex
+    to the polygon
+  * :meth:`~geographiclib.polygonarea.PolygonArea.AddEdge` add an edge
+    to the polygon
+  * :meth:`~geographiclib.polygonarea.PolygonArea.Compute` compute the
+    properties of the polygon
+  * :meth:`~geographiclib.polygonarea.PolygonArea.TestPoint` compute the
+    properties of the polygon with a tentative additional vertex
+  * :meth:`~geographiclib.polygonarea.PolygonArea.TestEdge` compute the
+    properties of the polygon with a tentative additional edge
+
+The public attributes for this class are
+
+  * :attr:`~geographiclib.polygonarea.PolygonArea.earth`
+    :attr:`~geographiclib.polygonarea.PolygonArea.polyline`
+    :attr:`~geographiclib.polygonarea.PolygonArea.area0`
+    :attr:`~geographiclib.polygonarea.PolygonArea.num`
+    :attr:`~geographiclib.polygonarea.PolygonArea.lat1`
+    :attr:`~geographiclib.polygonarea.PolygonArea.lon1`
+
+"""
+# polygonarea.py
+#
+# This is a rather literal translation of the GeographicLib::PolygonArea class
+# to python.  See the documentation for the C++ class for more information at
+#
+#    https://geographiclib.sourceforge.io/html/annotated.html
+#
+# The algorithms are derived in
+#
+#    Charles F. F. Karney,
+#    Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
+#    https://doi.org/10.1007/s00190-012-0578-z
+#    Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
+#
+# Copyright (c) Charles Karney (2011-2022) <charles@karney.com> and licensed
+# under the MIT/X11 License.  For more information, see
+# https://geographiclib.sourceforge.io/
+######################################################################
+
+import math
+from geographiclib.geomath import Math
+from geographiclib.accumulator import Accumulator
+from geographiclib.geodesic import Geodesic
+
+class PolygonArea:
+  """Area of a geodesic polygon"""
+
+  @staticmethod
+  def _transit(lon1, lon2):
+    """Count crossings of prime meridian for AddPoint."""
+    # Return 1 or -1 if crossing prime meridian in east or west direction.
+    # Otherwise return zero.
+    # Compute lon12 the same way as Geodesic::Inverse.
+    lon12, _ = Math.AngDiff(lon1, lon2)
+    lon1 = Math.AngNormalize(lon1)
+    lon2 = Math.AngNormalize(lon2)
+    return (1 if lon12 > 0 and ( lon1 < 0 <= lon2 or
+                                (lon1 > 0 and lon2 == 0))
+            else (-1 if lon12 < 0 and lon2 < 0 <= lon1 else 0))
+
+  @staticmethod
+  def _transitdirect(lon1, lon2):
+    """Count crossings of prime meridian for AddEdge."""
+    # We want to compute exactly
+    #   int(floor(lon2 / 360)) - int(floor(lon1 / 360))
+    lon1 = Math.remainder(lon1, 720.0); lon2 = Math.remainder(lon2, 720.0)
+    return ( (0 if 0 <= lon2 < 360 else 1) -
+             (0 if 0 <= lon1 < 360 else 1) )
+
+  @staticmethod
+  def _areareduceA(area, area0, crossings, reverse, sign):
+    """Reduce accumulator area to allowed range."""
+    area.Remainder(area0)
+    if crossings & 1:
+      area.Add( (1 if area.Sum() < 0 else -1) * area0/2 )
+    # area is with the clockwise sense.  If !reverse convert to
+    # counter-clockwise convention.
+    if not reverse: area.Negate()
+    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+    if sign:
+      if area.Sum() > area0/2:
+        area.Add( -area0 )
+      elif area.Sum() <= -area0/2:
+        area.Add(  area0 )
+    else:
+      if area.Sum() >= area0:
+        area.Add( -area0 )
+      elif area.Sum() < 0:
+        area.Add(  area0 )
+
+    return 0.0 + area.Sum()
+
+  @staticmethod
+  def _areareduceB(area, area0, crossings, reverse, sign):
+    """Reduce double area to allowed range."""
+    area = Math.remainder(area, area0)
+    if crossings & 1:
+      area += (1 if area < 0 else -1) * area0/2
+    # area is with the clockwise sense.  If !reverse convert to
+    # counter-clockwise convention.
+    if not reverse: area *= -1
+    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+    if sign:
+      if area > area0/2:
+        area -= area0
+      elif area <= -area0/2:
+        area += area0
+    else:
+      if area >= area0:
+        area -= area0
+      elif area < 0:
+        area += area0
+
+    return 0.0 + area
+
+  def __init__(self, earth, polyline = False):
+    """Construct a PolygonArea object
+
+    :param earth: a :class:`~geographiclib.geodesic.Geodesic` object
+    :param polyline: if true, treat object as a polyline instead of a polygon
+
+    Initially the polygon has no vertices.
+    """
+
+    self.earth = earth
+    """The geodesic object (readonly)"""
+    self.polyline = polyline
+    """Is this a polyline? (readonly)"""
+    self.area0 = 4 * math.pi * earth._c2
+    """The total area of the ellipsoid in meter^2 (readonly)"""
+    self._mask = (Geodesic.LATITUDE | Geodesic.LONGITUDE |
+                  Geodesic.DISTANCE |
+                  (Geodesic.EMPTY if self.polyline else
+                   Geodesic.AREA | Geodesic.LONG_UNROLL))
+    if not self.polyline: self._areasum = Accumulator()
+    self._perimetersum = Accumulator()
+    self.num = 0
+    """The current number of points in the polygon (readonly)"""
+    self.lat1 = math.nan
+    """The current latitude in degrees (readonly)"""
+    self.lon1 = math.nan
+    """The current longitude in degrees (readonly)"""
+    self._crossings = 0
+    self._lat0 = self._lon0 = math.nan
+
+  def Clear(self):
+    """Reset to empty polygon."""
+    self.num = 0
+    self._crossings = 0
+    if not self.polyline: self._areasum.Set(0)
+    self._perimetersum.Set(0)
+    self._lat0 = self._lon0 = self.lat1 = self.lon1 = math.nan
+
+  def AddPoint(self, lat, lon):
+    """Add the next vertex to the polygon
+
+    :param lat: the latitude of the point in degrees
+    :param lon: the longitude of the point in degrees
+
+    This adds an edge from the current vertex to the new vertex.
+    """
+
+    if self.num == 0:
+      self._lat0 = self.lat1 = lat
+      self._lon0 = self.lon1 = lon
+    else:
+      _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
+        self.lat1, self.lon1, lat, lon, self._mask)
+      self._perimetersum.Add(s12)
+      if not self.polyline:
+        self._areasum.Add(S12)
+        self._crossings += PolygonArea._transit(self.lon1, lon)
+      self.lat1 = lat
+      self.lon1 = lon
+    self.num += 1
+
+  def AddEdge(self, azi, s):
+    """Add the next edge to the polygon
+
+    :param azi: the azimuth at the current the point in degrees
+    :param s: the length of the edge in meters
+
+    This specifies the new vertex in terms of the edge from the current
+    vertex.
+
+    """
+
+    if self.num != 0:
+      _, lat, lon, _, _, _, _, _, S12 = self.earth._GenDirect(
+        self.lat1, self.lon1, azi, False, s, self._mask)
+      self._perimetersum.Add(s)
+      if not self.polyline:
+        self._areasum.Add(S12)
+        self._crossings += PolygonArea._transitdirect(self.lon1, lon)
+      self.lat1 = lat
+      self.lon1 = lon
+      self.num += 1
+
+  # return number, perimeter, area
+  def Compute(self, reverse = False, sign = True):
+    """Compute the properties of the polygon
+
+    :param reverse: if true then clockwise (instead of
+      counter-clockwise) traversal counts as a positive area
+    :param sign: if true then return a signed result for the area if the
+      polygon is traversed in the "wrong" direction instead of returning
+      the area for the rest of the earth
+    :return: a tuple of number, perimeter (meters), area (meters^2)
+
+    Arbitrarily complex polygons are allowed.  In the case of
+    self-intersecting polygons the area is accumulated "algebraically",
+    e.g., the areas of the 2 loops in a figure-8 polygon will partially
+    cancel.
+
+    If the object is a polygon (and not a polyline), the perimeter
+    includes the length of a final edge connecting the current point to
+    the initial point.  If the object is a polyline, then area is nan.
+
+    More points can be added to the polygon after this call.
+
+    """
+    if self.polyline: area = math.nan
+    if self.num < 2:
+      perimeter = 0.0
+      if not self.polyline: area = 0.0
+      return self.num, perimeter, area
+
+    if self.polyline:
+      perimeter = self._perimetersum.Sum()
+      return self.num, perimeter, area
+
+    _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
+      self.lat1, self.lon1, self._lat0, self._lon0, self._mask)
+    perimeter = self._perimetersum.Sum(s12)
+    tempsum = Accumulator(self._areasum)
+    tempsum.Add(S12)
+    crossings = self._crossings + PolygonArea._transit(self.lon1, self._lon0)
+    area = PolygonArea._areareduceA(tempsum, self.area0, crossings,
+                                    reverse, sign)
+    return self.num, perimeter, area
+
+  # return number, perimeter, area
+  def TestPoint(self, lat, lon, reverse = False, sign = True):
+    """Compute the properties for a tentative additional vertex
+
+    :param lat: the latitude of the point in degrees
+    :param lon: the longitude of the point in degrees
+    :param reverse: if true then clockwise (instead of
+      counter-clockwise) traversal counts as a positive area
+    :param sign: if true then return a signed result for the area if the
+      polygon is traversed in the "wrong" direction instead of returning
+      the area for the rest of the earth
+    :return: a tuple of number, perimeter (meters), area (meters^2)
+
+    """
+    if self.polyline: area = math.nan
+    if self.num == 0:
+      perimeter = 0.0
+      if not self.polyline: area = 0.0
+      return 1, perimeter, area
+
+    perimeter = self._perimetersum.Sum()
+    tempsum = 0.0 if self.polyline else self._areasum.Sum()
+    crossings = self._crossings; num = self.num + 1
+    for i in ([0] if self.polyline else [0, 1]):
+      _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
+        self.lat1 if i == 0 else lat, self.lon1 if i == 0 else lon,
+        self._lat0 if i != 0 else lat, self._lon0 if i != 0 else lon,
+        self._mask)
+      perimeter += s12
+      if not self.polyline:
+        tempsum += S12
+        crossings += PolygonArea._transit(self.lon1 if i == 0 else lon,
+                                         self._lon0 if i != 0 else lon)
+
+    if self.polyline:
+      return num, perimeter, area
+
+    area = PolygonArea._areareduceB(tempsum, self.area0, crossings,
+                                    reverse, sign)
+    return num, perimeter, area
+
+  # return num, perimeter, area
+  def TestEdge(self, azi, s, reverse = False, sign = True):
+    """Compute the properties for a tentative additional edge
+
+    :param azi: the azimuth at the current the point in degrees
+    :param s: the length of the edge in meters
+    :param reverse: if true then clockwise (instead of
+      counter-clockwise) traversal counts as a positive area
+    :param sign: if true then return a signed result for the area if the
+      polygon is traversed in the "wrong" direction instead of returning
+      the area for the rest of the earth
+    :return: a tuple of number, perimeter (meters), area (meters^2)
+
+    """
+
+    if self.num == 0:           # we don't have a starting point!
+      return 0, math.nan, math.nan
+    num = self.num + 1
+    perimeter = self._perimetersum.Sum() + s
+    if self.polyline:
+      return num, perimeter, math.nan
+
+    tempsum =  self._areasum.Sum()
+    crossings = self._crossings
+    _, lat, lon, _, _, _, _, _, S12 = self.earth._GenDirect(
+      self.lat1, self.lon1, azi, False, s, self._mask)
+    tempsum += S12
+    crossings += PolygonArea._transitdirect(self.lon1, lon)
+    _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
+      lat, lon, self._lat0, self._lon0, self._mask)
+    perimeter += s12
+    tempsum += S12
+    crossings += PolygonArea._transit(lon, self._lon0)
+
+    area = PolygonArea._areareduceB(tempsum, self.area0, crossings,
+                                    reverse, sign)
+    return num, perimeter, area

streamlit-venv/lib/python3.10/site-packages/geographiclib/test/__init__.py

@@ -0,0 +1,12 @@
+"""
+
+test_geodesic: test the geodesic routines from GeographicLib
+
+Run these tests with one of
+
+    python2 -m unittest -v geographiclib.test.test_geodesic
+    python3 -m unittest -v geographiclib.test.test_geodesic
+
+executed in this directory's parent directory.
+
+"""

streamlit-venv/lib/python3.10/site-packages/geographiclib/test/test_geodesic.py

@@ -0,0 +1,794 @@
+"""Geodesic tests"""
+
+import unittest
+import math
+
+from geographiclib.geodesic import Geodesic
+
+class GeodesicTest(unittest.TestCase):
+  """Geodesic test suite"""
+
+  testcases = [
+    [35.60777, -139.44815, 111.098748429560326,
+     -11.17491, -69.95921, 129.289270889708762,
+     8935244.5604818305, 80.50729714281974, 6273170.2055303837,
+     0.16606318447386067, 0.16479116945612937, 12841384694976.432],
+    [55.52454, 106.05087, 22.020059880982801,
+     77.03196, 197.18234, 109.112041110671519,
+     4105086.1713924406, 36.892740690445894, 3828869.3344387607,
+     0.80076349608092607, 0.80101006984201008, 61674961290615.615],
+    [-21.97856, 142.59065, -32.44456876433189,
+     41.84138, 98.56635, -41.84359951440466,
+     8394328.894657671, 75.62930491011522, 6161154.5773110616,
+     0.24816339233950381, 0.24930251203627892, -6637997720646.717],
+    [-66.99028, 112.2363, 173.73491240878403,
+     -12.70631, 285.90344, 2.512956620913668,
+     11150344.2312080241, 100.278634181155759, 6289939.5670446687,
+     -0.17199490274700385, -0.17722569526345708, -121287239862139.744],
+    [-17.42761, 173.34268, -159.033557661192928,
+     -15.84784, 5.93557, -20.787484651536988,
+     16076603.1631180673, 144.640108810286253, 3732902.1583877189,
+     -0.81273638700070476, -0.81299800519154474, 97825992354058.708],
+    [32.84994, 48.28919, 150.492927788121982,
+     -56.28556, 202.29132, 48.113449399816759,
+     16727068.9438164461, 150.565799985466607, 3147838.1910180939,
+     -0.87334918086923126, -0.86505036767110637, -72445258525585.010],
+    [6.96833, 52.74123, 92.581585386317712,
+     -7.39675, 206.17291, 90.721692165923907,
+     17102477.2496958388, 154.147366239113561, 2772035.6169917581,
+     -0.89991282520302447, -0.89986892177110739, -1311796973197.995],
+    [-50.56724, -16.30485, -105.439679907590164,
+     -33.56571, -94.97412, -47.348547835650331,
+     6455670.5118668696, 58.083719495371259, 5409150.7979815838,
+     0.53053508035997263, 0.52988722644436602, 41071447902810.047],
+    [-58.93002, -8.90775, 140.965397902500679,
+     -8.91104, 133.13503, 19.255429433416599,
+     11756066.0219864627, 105.755691241406877, 6151101.2270708536,
+     -0.26548622269867183, -0.27068483874510741, -86143460552774.735],
+    [-68.82867, -74.28391, 93.774347763114881,
+     -50.63005, -8.36685, 34.65564085411343,
+     3956936.926063544, 35.572254987389284, 3708890.9544062657,
+     0.81443963736383502, 0.81420859815358342, -41845309450093.787],
+    [-10.62672, -32.0898, -86.426713286747751,
+     5.883, -134.31681, -80.473780971034875,
+     11470869.3864563009, 103.387395634504061, 6184411.6622659713,
+     -0.23138683500430237, -0.23155097622286792, 4198803992123.548],
+    [-21.76221, 166.90563, 29.319421206936428,
+     48.72884, 213.97627, 43.508671946410168,
+     9098627.3986554915, 81.963476716121964, 6299240.9166992283,
+     0.13965943368590333, 0.14152969707656796, 10024709850277.476],
+    [-19.79938, -174.47484, 71.167275780171533,
+     -11.99349, -154.35109, 65.589099775199228,
+     2319004.8601169389, 20.896611684802389, 2267960.8703918325,
+     0.93427001867125849, 0.93424887135032789, -3935477535005.785],
+    [-11.95887, -116.94513, 92.712619830452549,
+     4.57352, 7.16501, 78.64960934409585,
+     13834722.5801401374, 124.688684161089762, 5228093.177931598,
+     -0.56879356755666463, -0.56918731952397221, -9919582785894.853],
+    [-87.85331, 85.66836, -65.120313040242748,
+     66.48646, 16.09921, -4.888658719272296,
+     17286615.3147144645, 155.58592449699137, 2635887.4729110181,
+     -0.90697975771398578, -0.91095608883042767, 42667211366919.534],
+    [1.74708, 128.32011, -101.584843631173858,
+     -11.16617, 11.87109, -86.325793296437476,
+     12942901.1241347408, 116.650512484301857, 5682744.8413270572,
+     -0.44857868222697644, -0.44824490340007729, 10763055294345.653],
+    [-25.72959, -144.90758, -153.647468693117198,
+     -57.70581, -269.17879, -48.343983158876487,
+     9413446.7452453107, 84.664533838404295, 6356176.6898881281,
+     0.09492245755254703, 0.09737058264766572, 74515122850712.444],
+    [-41.22777, 122.32875, 14.285113402275739,
+     -7.57291, 130.37946, 10.805303085187369,
+     3812686.035106021, 34.34330804743883, 3588703.8812128856,
+     0.82605222593217889, 0.82572158200920196, -2456961531057.857],
+    [11.01307, 138.25278, 79.43682622782374,
+     6.62726, 247.05981, 103.708090215522657,
+     11911190.819018408, 107.341669954114577, 6070904.722786735,
+     -0.29767608923657404, -0.29785143390252321, 17121631423099.696],
+    [-29.47124, 95.14681, -163.779130441688382,
+     -27.46601, -69.15955, -15.909335945554969,
+     13487015.8381145492, 121.294026715742277, 5481428.9945736388,
+     -0.51527225545373252, -0.51556587964721788, 104679964020340.318]]
+
+  def test_inverse(self):
+    """Helper function for testing inverse calculation"""
+    for l in GeodesicTest.testcases:
+      (lat1, lon1, azi1, lat2, lon2, azi2,
+       s12, a12, m12, M12, M21, S12) = l
+      inv = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2,
+                     Geodesic.ALL | Geodesic.LONG_UNROLL)
+      self.assertAlmostEqual(lon2, inv["lon2"], delta = 1e-13)
+      self.assertAlmostEqual(azi1, inv["azi1"], delta = 1e-13)
+      self.assertAlmostEqual(azi2, inv["azi2"], delta = 1e-13)
+      self.assertAlmostEqual(s12, inv["s12"], delta = 1e-8)
+      self.assertAlmostEqual(a12, inv["a12"], delta = 1e-13)
+      self.assertAlmostEqual(m12, inv["m12"], delta = 1e-8)
+      self.assertAlmostEqual(M12, inv["M12"], delta = 1e-15)
+      self.assertAlmostEqual(M21, inv["M21"], delta = 1e-15)
+      self.assertAlmostEqual(S12, inv["S12"], delta = 0.1)
+
+  def test_direct(self):
+    """Helper function for testing direct calculation"""
+    for l in GeodesicTest.testcases:
+      (lat1, lon1, azi1, lat2, lon2, azi2,
+       s12, a12, m12, M12, M21, S12) = l
+      direct = Geodesic.WGS84.Direct(lat1, lon1, azi1, s12,
+                                     Geodesic.ALL | Geodesic.LONG_UNROLL)
+      self.assertAlmostEqual(lat2, direct["lat2"], delta = 1e-13)
+      self.assertAlmostEqual(lon2, direct["lon2"], delta = 1e-13)
+      self.assertAlmostEqual(azi2, direct["azi2"], delta = 1e-13)
+      self.assertAlmostEqual(a12, direct["a12"], delta = 1e-13)
+      self.assertAlmostEqual(m12, direct["m12"], delta = 1e-8)
+      self.assertAlmostEqual(M12, direct["M12"], delta = 1e-15)
+      self.assertAlmostEqual(M21, direct["M21"], delta = 1e-15)
+      self.assertAlmostEqual(S12, direct["S12"], delta = 0.1)
+
+  def test_arcdirect(self):
+    """Helper function for testing direct calculation with arc length"""
+    for l in GeodesicTest.testcases:
+      (lat1, lon1, azi1, lat2, lon2, azi2,
+       s12, a12, m12, M12, M21, S12) = l
+      direct = Geodesic.WGS84.ArcDirect(lat1, lon1, azi1, a12,
+                                        Geodesic.ALL | Geodesic.LONG_UNROLL)
+      self.assertAlmostEqual(lat2, direct["lat2"], delta = 1e-13)
+      self.assertAlmostEqual(lon2, direct["lon2"], delta = 1e-13)
+      self.assertAlmostEqual(azi2, direct["azi2"], delta = 1e-13)
+      self.assertAlmostEqual(s12, direct["s12"], delta = 1e-8)
+      self.assertAlmostEqual(m12, direct["m12"], delta = 1e-8)
+      self.assertAlmostEqual(M12, direct["M12"], delta = 1e-15)
+      self.assertAlmostEqual(M21, direct["M21"], delta = 1e-15)
+      self.assertAlmostEqual(S12, direct["S12"], delta = 0.1)
+
+class GeodSolveTest(unittest.TestCase):
+  """GeodSolve tests"""
+
+  def test_GeodSolve0(self):
+    """GeodSolve0"""
+    inv = Geodesic.WGS84.Inverse(40.6, -73.8, 49.01666667, 2.55)
+    self.assertAlmostEqual(inv["azi1"], 53.47022, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 111.59367, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 5853226, delta = 0.5)
+
+  def test_GeodSolve1(self):
+    """GeodSolve1"""
+    direct = Geodesic.WGS84.Direct(40.63972222, -73.77888889, 53.5, 5850e3)
+    self.assertAlmostEqual(direct["lat2"], 49.01467, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], 2.56106, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["azi2"], 111.62947, delta = 0.5e-5)
+
+  def test_GeodSolve2(self):
+    """Check fix for antipodal prolate bug found 2010-09-04"""
+    geod = Geodesic(6.4e6, -1/150.0)
+    inv = geod.Inverse(0.07476, 0, -0.07476, 180)
+    self.assertAlmostEqual(inv["azi1"], 90.00078, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 90.00078, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20106193, delta = 0.5)
+    inv = geod.Inverse(0.1, 0, -0.1, 180)
+    self.assertAlmostEqual(inv["azi1"], 90.00105, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 90.00105, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20106193, delta = 0.5)
+
+  def test_GeodSolve4(self):
+    """Check fix for short line bug found 2010-05-21"""
+    inv = Geodesic.WGS84.Inverse(36.493349428792, 0,
+                                 36.49334942879201, 0.0000008)
+    self.assertAlmostEqual(inv["s12"], 0.072, delta = 0.5e-3)
+
+  def test_GeodSolve5(self):
+    """Check fix for point2=pole bug found 2010-05-03"""
+    direct = Geodesic.WGS84.Direct(0.01777745589997, 30, 0, 10e6)
+    self.assertAlmostEqual(direct["lat2"], 90, delta = 0.5e-5)
+    if direct["lon2"] < 0:
+      self.assertAlmostEqual(direct["lon2"], -150, delta = 0.5e-5)
+      self.assertAlmostEqual(abs(direct["azi2"]), 180, delta = 0.5e-5)
+    else:
+      self.assertAlmostEqual(direct["lon2"], 30, delta = 0.5e-5)
+      self.assertAlmostEqual(direct["azi2"], 0, delta = 0.5e-5)
+
+  def test_GeodSolve6(self):
+    """Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4
+    x86 -O3).  Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1)."""
+    inv = Geodesic.WGS84.Inverse(88.202499451857, 0,
+                                 -88.202499451857, 179.981022032992859592)
+    self.assertAlmostEqual(inv["s12"], 20003898.214, delta = 0.5e-3)
+    inv = Geodesic.WGS84.Inverse(89.262080389218, 0,
+                                 -89.262080389218, 179.992207982775375662)
+    self.assertAlmostEqual(inv["s12"], 20003925.854, delta = 0.5e-3)
+    inv = Geodesic.WGS84.Inverse(89.333123580033, 0,
+                                 -89.333123580032997687, 179.99295812360148422)
+    self.assertAlmostEqual(inv["s12"], 20003926.881, delta = 0.5e-3)
+
+  def test_GeodSolve9(self):
+    """Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3)"""
+    inv = Geodesic.WGS84.Inverse(56.320923501171, 0,
+                                 -56.320923501171, 179.664747671772880215)
+    self.assertAlmostEqual(inv["s12"], 19993558.287, delta = 0.5e-3)
+
+  def test_GeodSolve10(self):
+    """Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio
+    10 rel + debug)"""
+    inv = Geodesic.WGS84.Inverse(52.784459512564, 0,
+                                 -52.784459512563990912, 179.634407464943777557)
+    self.assertAlmostEqual(inv["s12"], 19991596.095, delta = 0.5e-3)
+
+  def test_GeodSolve11(self):
+    """Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio
+    10 rel + debug)"""
+    inv = Geodesic.WGS84.Inverse(48.522876735459, 0,
+                                 -48.52287673545898293, 179.599720456223079643)
+    self.assertAlmostEqual(inv["s12"], 19989144.774, delta = 0.5e-3)
+
+  def test_GeodSolve12(self):
+    """Check fix for inverse geodesics on extreme prolate/oblate
+    ellipsoids Reported 2012-08-29 Stefan Guenther
+    <stefan.gunther@embl.de>; fixed 2012-10-07"""
+    geod = Geodesic(89.8, -1.83)
+    inv = geod.Inverse(0, 0, -10, 160)
+    self.assertAlmostEqual(inv["azi1"], 120.27, delta = 1e-2)
+    self.assertAlmostEqual(inv["azi2"], 105.15, delta = 1e-2)
+    self.assertAlmostEqual(inv["s12"], 266.7, delta = 1e-1)
+
+  def test_GeodSolve14(self):
+    """Check fix for inverse ignoring lon12 = nan"""
+    inv = Geodesic.WGS84.Inverse(0, 0, 1, math.nan)
+    self.assertTrue(math.isnan(inv["azi1"]))
+    self.assertTrue(math.isnan(inv["azi2"]))
+    self.assertTrue(math.isnan(inv["s12"]))
+
+  def test_GeodSolve15(self):
+    """Initial implementation of Math::eatanhe was wrong for e^2 < 0.  This
+    checks that this is fixed."""
+    geod = Geodesic(6.4e6, -1/150.0)
+    direct = geod.Direct(1, 2, 3, 4, Geodesic.AREA)
+    self.assertAlmostEqual(direct["S12"], 23700, delta = 0.5)
+
+  def test_GeodSolve17(self):
+    """Check fix for LONG_UNROLL bug found on 2015-05-07"""
+    direct = Geodesic.WGS84.Direct(40, -75, -10, 2e7,
+                                   Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], -39, delta = 1)
+    self.assertAlmostEqual(direct["lon2"], -254, delta = 1)
+    self.assertAlmostEqual(direct["azi2"], -170, delta = 1)
+    line = Geodesic.WGS84.Line(40, -75, -10)
+    direct = line.Position(2e7, Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], -39, delta = 1)
+    self.assertAlmostEqual(direct["lon2"], -254, delta = 1)
+    self.assertAlmostEqual(direct["azi2"], -170, delta = 1)
+    direct = Geodesic.WGS84.Direct(40, -75, -10, 2e7)
+    self.assertAlmostEqual(direct["lat2"], -39, delta = 1)
+    self.assertAlmostEqual(direct["lon2"], 105, delta = 1)
+    self.assertAlmostEqual(direct["azi2"], -170, delta = 1)
+    direct = line.Position(2e7)
+    self.assertAlmostEqual(direct["lat2"], -39, delta = 1)
+    self.assertAlmostEqual(direct["lon2"], 105, delta = 1)
+    self.assertAlmostEqual(direct["azi2"], -170, delta = 1)
+
+  def test_GeodSolve26(self):
+    """Check 0/0 problem with area calculation on sphere 2015-09-08"""
+    geod = Geodesic(6.4e6, 0)
+    inv = geod.Inverse(1, 2, 3, 4, Geodesic.AREA)
+    self.assertAlmostEqual(inv["S12"], 49911046115.0, delta = 0.5)
+
+  def test_GeodSolve28(self):
+    """Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in
+    Java implementation fixed on 2015-05-19)."""
+    geod = Geodesic(6.4e6, 0.1)
+    direct = geod.Direct(1, 2, 10, 5e6)
+    self.assertAlmostEqual(direct["a12"], 48.55570690, delta = 0.5e-8)
+
+  def test_GeodSolve29(self):
+    """Check longitude unrolling with inverse calculation 2015-09-16"""
+    direct = Geodesic.WGS84.Inverse(0, 539, 0, 181)
+    self.assertAlmostEqual(direct["lon1"], 179, delta = 1e-10)
+    self.assertAlmostEqual(direct["lon2"], -179, delta = 1e-10)
+    self.assertAlmostEqual(direct["s12"], 222639, delta = 0.5)
+    direct = Geodesic.WGS84.Inverse(0, 539, 0, 181,
+                                    Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lon1"], 539, delta = 1e-10)
+    self.assertAlmostEqual(direct["lon2"], 541, delta = 1e-10)
+    self.assertAlmostEqual(direct["s12"], 222639, delta = 0.5)
+
+  def test_GeodSolve33(self):
+    """Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in
+    Octave -- sind(-0.0) = +0.0 -- and in some version of Visual
+    Studio -- fmod(-0.0, 360.0) = +0.0."""
+    inv = Geodesic.WGS84.Inverse(0, 0, 0, 179)
+    self.assertAlmostEqual(inv["azi1"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 19926189, delta = 0.5)
+    inv = Geodesic.WGS84.Inverse(0, 0, 0, 179.5)
+    self.assertAlmostEqual(inv["azi1"], 55.96650, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 124.03350, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 19980862, delta = 0.5)
+    inv = Geodesic.WGS84.Inverse(0, 0, 0, 180)
+    self.assertAlmostEqual(inv["azi1"], 0.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(inv["azi2"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20003931, delta = 0.5)
+    inv = Geodesic.WGS84.Inverse(0, 0, 1, 180)
+    self.assertAlmostEqual(inv["azi1"], 0.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(inv["azi2"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 19893357, delta = 0.5)
+    geod = Geodesic(6.4e6, 0)
+    inv = geod.Inverse(0, 0, 0, 179)
+    self.assertAlmostEqual(inv["azi1"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 19994492, delta = 0.5)
+    inv = geod.Inverse(0, 0, 0, 180)
+    self.assertAlmostEqual(inv["azi1"], 0.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(inv["azi2"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20106193, delta = 0.5)
+    inv = geod.Inverse(0, 0, 1, 180)
+    self.assertAlmostEqual(inv["azi1"], 0.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(inv["azi2"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 19994492, delta = 0.5)
+    geod = Geodesic(6.4e6, -1/300.0)
+    inv = geod.Inverse(0, 0, 0, 179)
+    self.assertAlmostEqual(inv["azi1"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 19994492, delta = 0.5)
+    inv = geod.Inverse(0, 0, 0, 180)
+    self.assertAlmostEqual(inv["azi1"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 90.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20106193, delta = 0.5)
+    inv = geod.Inverse(0, 0, 0.5, 180)
+    self.assertAlmostEqual(inv["azi1"], 33.02493, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["azi2"], 146.97364, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20082617, delta = 0.5)
+    inv = geod.Inverse(0, 0, 1, 180)
+    self.assertAlmostEqual(inv["azi1"], 0.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(inv["azi2"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(inv["s12"], 20027270, delta = 0.5)
+
+  def test_GeodSolve55(self):
+    """Check fix for nan + point on equator or pole not returning all nans in
+    Geodesic::Inverse, found 2015-09-23."""
+    inv = Geodesic.WGS84.Inverse(math.nan, 0, 0, 90)
+    self.assertTrue(math.isnan(inv["azi1"]))
+    self.assertTrue(math.isnan(inv["azi2"]))
+    self.assertTrue(math.isnan(inv["s12"]))
+    inv = Geodesic.WGS84.Inverse(math.nan, 0, 90, 9)
+    self.assertTrue(math.isnan(inv["azi1"]))
+    self.assertTrue(math.isnan(inv["azi2"]))
+    self.assertTrue(math.isnan(inv["s12"]))
+
+  def test_GeodSolve59(self):
+    """Check for points close with longitudes close to 180 deg apart."""
+    inv = Geodesic.WGS84.Inverse(5, 0.00000000000001, 10, 180)
+    self.assertAlmostEqual(inv["azi1"], 0.000000000000035, delta = 1.5e-14)
+    self.assertAlmostEqual(inv["azi2"], 179.99999999999996, delta = 1.5e-14)
+    self.assertAlmostEqual(inv["s12"], 18345191.174332713, delta = 5e-9)
+
+  def test_GeodSolve61(self):
+    """Make sure small negative azimuths are west-going"""
+    direct = Geodesic.WGS84.Direct(45, 0, -0.000000000000000003, 1e7,
+                                   Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], 45.30632, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -180, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(direct["azi2"]), 180, delta = 0.5e-5)
+    line = Geodesic.WGS84.InverseLine(45, 0, 80, -0.000000000000000003)
+    direct = line.Position(1e7, Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], 45.30632, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -180, delta = 0.5e-5)
+    self.assertAlmostEqual(abs(direct["azi2"]), 180, delta = 0.5e-5)
+
+  def test_GeodSolve65(self):
+    """Check for bug in east-going check in GeodesicLine (needed to check for
+    sign of 0) and sign error in area calculation due to a bogus override
+    of the code for alp12.  Found/fixed on 2015-12-19."""
+    line = Geodesic.WGS84.InverseLine(30, -0.000000000000000001, -31, 180,
+                                      Geodesic.ALL)
+    direct = line.Position(1e7, Geodesic.ALL | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat1"], 30.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon1"], -0.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(abs(direct["azi1"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lat2"], -60.23169 , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -0.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(abs(direct["azi2"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["s12"] , 10000000  , delta = 0.5)
+    self.assertAlmostEqual(direct["a12"] , 90.06544  , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["m12"] , 6363636   , delta = 0.5)
+    self.assertAlmostEqual(direct["M12"] , -0.0012834, delta = 0.5e7)
+    self.assertAlmostEqual(direct["M21"] , 0.0013749 , delta = 0.5e-7)
+    self.assertAlmostEqual(direct["S12"] , 0         , delta = 0.5)
+    direct = line.Position(2e7, Geodesic.ALL | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat1"], 30.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon1"], -0.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(abs(direct["azi1"]), 180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lat2"], -30.03547 , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["azi2"], -0.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["s12"] , 20000000  , delta = 0.5)
+    self.assertAlmostEqual(direct["a12"] , 179.96459 , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["m12"] , 54342     , delta = 0.5)
+    self.assertAlmostEqual(direct["M12"] , -1.0045592, delta = 0.5e7)
+    self.assertAlmostEqual(direct["M21"] , -0.9954339, delta = 0.5e-7)
+    self.assertAlmostEqual(direct["S12"] , 127516405431022.0, delta = 0.5)
+
+  def test_GeodSolve66(self):
+    """Check for InverseLine if line is slightly west of S and that s13 is
+    correctly set."""
+    line = Geodesic.WGS84.InverseLine(-5, -0.000000000000002, -10, 180)
+    direct = line.Position(2e7, Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], 4.96445   , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -180.00000, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["azi2"], -0.00000  , delta = 0.5e-5)
+    direct = line.Position(0.5 * line.s13,
+                           Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], -87.52461 , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -0.00000  , delta = 0.5e-5)
+    self.assertAlmostEqual(direct["azi2"], -180.00000, delta = 0.5e-5)
+
+  def test_GeodSolve71(self):
+    """Check that DirectLine sets s13."""
+    line = Geodesic.WGS84.DirectLine(1, 2, 45, 1e7)
+    direct = line.Position(0.5 * line.s13,
+                           Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+    self.assertAlmostEqual(direct["lat2"], 30.92625, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], 37.54640, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["azi2"], 55.43104, delta = 0.5e-5)
+
+  def test_GeodSolve73(self):
+    """Check for backwards from the pole bug reported by Anon on 2016-02-13.
+    This only affected the Java implementation.  It was introduced in Java
+    version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17.
+    Also the + sign on azi2 is a check on the normalizing of azimuths
+    (converting -0.0 to +0.0)."""
+    direct = Geodesic.WGS84.Direct(90, 10, 180, -1e6)
+    self.assertAlmostEqual(direct["lat2"], 81.04623, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["lon2"], -170, delta = 0.5e-5)
+    self.assertAlmostEqual(direct["azi2"], 0, delta = 0.5e-5)
+    self.assertTrue(math.copysign(1, direct["azi2"]) > 0)
+
+  def test_GeodSolve74(self):
+    """Check fix for inaccurate areas, bug introduced in v1.46, fixed
+    2015-10-16."""
+    inv = Geodesic.WGS84.Inverse(54.1589, 15.3872, 54.1591, 15.3877,
+                                 Geodesic.ALL)
+    self.assertAlmostEqual(inv["azi1"], 55.723110355, delta = 5e-9)
+    self.assertAlmostEqual(inv["azi2"], 55.723515675, delta = 5e-9)
+    self.assertAlmostEqual(inv["s12"],  39.527686385, delta = 5e-9)
+    self.assertAlmostEqual(inv["a12"],   0.000355495, delta = 5e-9)
+    self.assertAlmostEqual(inv["m12"],  39.527686385, delta = 5e-9)
+    self.assertAlmostEqual(inv["M12"],   0.999999995, delta = 5e-9)
+    self.assertAlmostEqual(inv["M21"],   0.999999995, delta = 5e-9)
+    self.assertAlmostEqual(inv["S12"], 286698586.30197, delta = 5e-4)
+
+  def test_GeodSolve76(self):
+    """The distance from Wellington and Salamanca (a classic failure of
+    Vincenty)"""
+    inv = Geodesic.WGS84.Inverse(-(41+19/60.0), 174+49/60.0,
+                                 40+58/60.0, -(5+30/60.0))
+    self.assertAlmostEqual(inv["azi1"], 160.39137649664, delta = 0.5e-11)
+    self.assertAlmostEqual(inv["azi2"],  19.50042925176, delta = 0.5e-11)
+    self.assertAlmostEqual(inv["s12"],  19960543.857179, delta = 0.5e-6)
+
+  def test_GeodSolve78(self):
+    """An example where the NGS calculator fails to converge"""
+    inv = Geodesic.WGS84.Inverse(27.2, 0.0, -27.1, 179.5)
+    self.assertAlmostEqual(inv["azi1"],  45.82468716758, delta = 0.5e-11)
+    self.assertAlmostEqual(inv["azi2"], 134.22776532670, delta = 0.5e-11)
+    self.assertAlmostEqual(inv["s12"],  19974354.765767, delta = 0.5e-6)
+
+  def test_GeodSolve80(self):
+    """Some tests to add code coverage: computing scale in special cases + zero
+    length geodesic (includes GeodSolve80 - GeodSolve83) + using an incapable
+    line."""
+    inv = Geodesic.WGS84.Inverse(0, 0, 0, 90, Geodesic.GEODESICSCALE)
+    self.assertAlmostEqual(inv["M12"], -0.00528427534, delta = 0.5e-10)
+    self.assertAlmostEqual(inv["M21"], -0.00528427534, delta = 0.5e-10)
+
+    inv = Geodesic.WGS84.Inverse(0, 0, 1e-6, 1e-6, Geodesic.GEODESICSCALE)
+    self.assertAlmostEqual(inv["M12"], 1, delta = 0.5e-10)
+    self.assertAlmostEqual(inv["M21"], 1, delta = 0.5e-10)
+
+    inv = Geodesic.WGS84.Inverse(20.001, 0, 20.001, 0, Geodesic.ALL)
+    self.assertAlmostEqual(inv["a12"], 0, delta = 1e-13)
+    self.assertAlmostEqual(inv["s12"], 0, delta = 1e-8)
+    self.assertAlmostEqual(inv["azi1"], 180, delta = 1e-13)
+    self.assertAlmostEqual(inv["azi2"], 180, delta = 1e-13)
+    self.assertAlmostEqual(inv["m12"], 0, delta =  1e-8)
+    self.assertAlmostEqual(inv["M12"], 1, delta = 1e-15)
+    self.assertAlmostEqual(inv["M21"], 1, delta = 1e-15)
+    self.assertAlmostEqual(inv["S12"], 0, delta = 1e-10)
+    self.assertTrue(math.copysign(1, inv["a12"]) > 0)
+    self.assertTrue(math.copysign(1, inv["s12"]) > 0)
+    self.assertTrue(math.copysign(1, inv["m12"]) > 0)
+
+    inv = Geodesic.WGS84.Inverse(90, 0, 90, 180, Geodesic.ALL)
+    self.assertAlmostEqual(inv["a12"], 0, delta = 1e-13)
+    self.assertAlmostEqual(inv["s12"], 0, delta = 1e-8)
+    self.assertAlmostEqual(inv["azi1"], 0, delta = 1e-13)
+    self.assertAlmostEqual(inv["azi2"], 180, delta = 1e-13)
+    self.assertAlmostEqual(inv["m12"], 0, delta = 1e-8)
+    self.assertAlmostEqual(inv["M12"], 1, delta = 1e-15)
+    self.assertAlmostEqual(inv["M21"], 1, delta = 1e-15)
+    self.assertAlmostEqual(inv["S12"], 127516405431022.0, delta = 0.5)
+
+    # An incapable line which can't take distance as input
+    line = Geodesic.WGS84.Line(1, 2, 90, Geodesic.LATITUDE)
+    direct = line.Position(1000, Geodesic.EMPTY)
+    self.assertTrue(math.isnan(direct["a12"]))
+
+  def test_GeodSolve84(self):
+    """Tests for python implementation to check fix for range errors with
+    {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve91)."""
+    direct = Geodesic.WGS84.Direct(0, 0, 90, math.inf)
+    self.assertTrue(math.isnan(direct["lat2"]))
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(math.isnan(direct["azi2"]))
+    direct = Geodesic.WGS84.Direct(0, 0, 90, math.nan)
+    self.assertTrue(math.isnan(direct["lat2"]))
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(math.isnan(direct["azi2"]))
+    direct = Geodesic.WGS84.Direct(0, 0, math.inf, 1000)
+    self.assertTrue(math.isnan(direct["lat2"]))
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(math.isnan(direct["azi2"]))
+    direct = Geodesic.WGS84.Direct(0, 0, math.nan, 1000)
+    self.assertTrue(math.isnan(direct["lat2"]))
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(math.isnan(direct["azi2"]))
+    direct = Geodesic.WGS84.Direct(0, math.inf, 90, 1000)
+    self.assertTrue(direct["lat1"] == 0)
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(direct["azi2"] == 90)
+    direct = Geodesic.WGS84.Direct(0, math.nan, 90, 1000)
+    self.assertTrue(direct["lat1"] == 0)
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(direct["azi2"] == 90)
+    direct = Geodesic.WGS84.Direct(math.inf, 0, 90, 1000)
+    self.assertTrue(math.isnan(direct["lat2"]))
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(math.isnan(direct["azi2"]))
+    direct = Geodesic.WGS84.Direct(math.nan, 0, 90, 1000)
+    self.assertTrue(math.isnan(direct["lat2"]))
+    self.assertTrue(math.isnan(direct["lon2"]))
+    self.assertTrue(math.isnan(direct["azi2"]))
+
+  def test_GeodSolve92(self):
+    """Check fix for inaccurate hypot with python 3.[89].  Problem reported
+    by agdhruv https://github.com/geopy/geopy/issues/466 ; see
+    https://bugs.python.org/issue43088"""
+    inv = Geodesic.WGS84.Inverse(37.757540000000006, -122.47018,
+                                 37.75754,           -122.470177)
+    self.assertAlmostEqual(inv["azi1"], 89.99999923, delta = 1e-7  )
+    self.assertAlmostEqual(inv["azi2"], 90.00000106, delta = 1e-7  )
+    self.assertAlmostEqual(inv["s12"],   0.264,      delta = 0.5e-3)
+
+  def test_GeodSolve94(self):
+    """Check fix for lat2 = nan being treated as lat2 = 0 (bug found
+       2021-07-26)"""
+    inv = Geodesic.WGS84.Inverse(0, 0, math.nan, 90)
+    self.assertTrue(math.isnan(inv["azi1"]))
+    self.assertTrue(math.isnan(inv["azi2"]))
+    self.assertTrue(math.isnan(inv["s12"]))
+
+  def test_GeodSolve96(self):
+    """Failure with long doubles found with test case from Nowak + Nowak Da
+       Costa (2022).  Problem was using somg12 > 1 as a test that it needed
+       to be set when roundoff could result in somg12 slightly bigger that 1.
+       Found + fixed 2022-03-30."""
+    geod = Geodesic(6378137, 1/298.257222101)
+    inv = geod.Inverse(0, 0, 60.0832522871723, 89.8492185074635, Geodesic.AREA)
+    self.assertAlmostEqual(inv["S12"], 42426932221845, delta = 0.5)
+
+class PlanimeterTest(unittest.TestCase):
+  """Planimeter tests"""
+
+  polygon = Geodesic.WGS84.Polygon(False)
+  polyline = Geodesic.WGS84.Polygon(True)
+
+  @staticmethod
+  def Planimeter(points):
+    """Helper function for polygons"""
+    PlanimeterTest.polygon.Clear()
+    for p in points:
+      PlanimeterTest.polygon.AddPoint(p[0], p[1])
+    return PlanimeterTest.polygon.Compute(False, True)
+
+  @staticmethod
+  def PolyLength(points):
+    """Helper function for polylines"""
+    PlanimeterTest.polyline.Clear()
+    for p in points:
+      PlanimeterTest.polyline.AddPoint(p[0], p[1])
+    return PlanimeterTest.polyline.Compute(False, True)
+
+  def test_Planimeter0(self):
+    """Check fix for pole-encircling bug found 2011-03-16"""
+    points = [[89, 0], [89, 90], [89, 180], [89, 270]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 631819.8745, delta = 1e-4)
+    self.assertAlmostEqual(area, 24952305678.0, delta = 1)
+    points = [[-89, 0], [-89, 90], [-89, 180], [-89, 270]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 631819.8745, delta = 1e-4)
+    self.assertAlmostEqual(area, -24952305678.0, delta = 1)
+
+    points = [[0, -1], [-1, 0], [0, 1], [1, 0]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 627598.2731, delta = 1e-4)
+    self.assertAlmostEqual(area, 24619419146.0, delta = 1)
+
+    points = [[90, 0], [0, 0], [0, 90]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 30022685, delta = 1)
+    self.assertAlmostEqual(area, 63758202715511.0, delta = 1)
+    _, perimeter, area = PlanimeterTest.PolyLength(points)
+    self.assertAlmostEqual(perimeter, 20020719, delta = 1)
+    self.assertTrue(math.isnan(area))
+
+  def test_Planimeter5(self):
+    """Check fix for Planimeter pole crossing bug found 2011-06-24"""
+    points = [[89, 0.1], [89, 90.1], [89, -179.9]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 539297, delta = 1)
+    self.assertAlmostEqual(area, 12476152838.5, delta = 1)
+
+  def test_Planimeter6(self):
+    """Check fix for Planimeter lon12 rounding bug found 2012-12-03"""
+    points = [[9, -0.00000000000001], [9, 180], [9, 0]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 36026861, delta = 1)
+    self.assertAlmostEqual(area, 0, delta = 1)
+    points = [[9, 0.00000000000001], [9, 0], [9, 180]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 36026861, delta = 1)
+    self.assertAlmostEqual(area, 0, delta = 1)
+    points = [[9, 0.00000000000001], [9, 180], [9, 0]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 36026861, delta = 1)
+    self.assertAlmostEqual(area, 0, delta = 1)
+    points = [[9, -0.00000000000001], [9, 0], [9, 180]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 36026861, delta = 1)
+    self.assertAlmostEqual(area, 0, delta = 1)
+
+  def test_Planimeter12(self):
+    """Area of arctic circle (not really -- adjunct to rhumb-area test)"""
+    points = [[66.562222222, 0], [66.562222222, 180], [66.562222222, 360]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 10465729, delta = 1)
+    self.assertAlmostEqual(area, 0, delta = 1)
+
+  def test_Planimeter12r(self):
+    """Reverse area of arctic circle"""
+    points = [[66.562222222, -0], [66.562222222, -180], [66.562222222, -360]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 10465729, delta = 1)
+    self.assertAlmostEqual(area, 0, delta = 1)
+
+  def test_Planimeter13(self):
+    """Check encircling pole twice"""
+    points = [[89,-360], [89,-240], [89,-120], [89,0], [89,120], [89,240]]
+    _, perimeter, area = PlanimeterTest.Planimeter(points)
+    self.assertAlmostEqual(perimeter, 1160741, delta = 1)
+    self.assertAlmostEqual(area, 32415230256.0, delta = 1)
+
+  def test_Planimeter15(self):
+    """Coverage tests, includes Planimeter15 - Planimeter18 (combinations of
+    reverse and sign) + calls to testpoint, testedge."""
+    lat = [2, 1, 3]
+    lon = [1, 2, 3]
+    r = 18454562325.45119
+    a0 = 510065621724088.5093   # ellipsoid area
+    PlanimeterTest.polygon.Clear()
+    PlanimeterTest.polygon.AddPoint(lat[0], lon[0])
+    PlanimeterTest.polygon.AddPoint(lat[1], lon[1])
+    _, _, area = PlanimeterTest.polygon.TestPoint(lat[2], lon[2], False, True)
+    self.assertAlmostEqual(area, r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.TestPoint(lat[2], lon[2], False, False)
+    self.assertAlmostEqual(area, r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.TestPoint(lat[2], lon[2], True, True)
+    self.assertAlmostEqual(area, -r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.TestPoint(lat[2], lon[2], True, False)
+    self.assertAlmostEqual(area, a0-r, delta = 0.5)
+    inv = Geodesic.WGS84.Inverse(lat[1], lon[1], lat[2], lon[2])
+    azi1 = inv["azi1"]
+    s12 = inv["s12"]
+    _, _, area = PlanimeterTest.polygon.TestEdge(azi1, s12, False, True)
+    self.assertAlmostEqual(area, r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.TestEdge(azi1, s12, False, False)
+    self.assertAlmostEqual(area, r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.TestEdge(azi1, s12, True, True)
+    self.assertAlmostEqual(area, -r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.TestEdge(azi1, s12, True, False)
+    self.assertAlmostEqual(area, a0-r, delta = 0.5)
+    PlanimeterTest.polygon.AddPoint(lat[2], lon[2])
+    _, _, area = PlanimeterTest.polygon.Compute(False, True)
+    self.assertAlmostEqual(area, r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.Compute(False, False)
+    self.assertAlmostEqual(area, r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.Compute(True, True)
+    self.assertAlmostEqual(area, -r, delta = 0.5)
+    _, _, area = PlanimeterTest.polygon.Compute(True, False)
+    self.assertAlmostEqual(area, a0-r, delta = 0.5)
+
+  def test_Planimeter19(self):
+    """Coverage tests, includes Planimeter19 - Planimeter20 (degenerate
+    polygons) + extra cases."""
+    PlanimeterTest.polygon.Clear()
+    _, perimeter, area = PlanimeterTest.polygon.Compute(False, True)
+    self.assertTrue(area == 0)
+    self.assertTrue(perimeter == 0)
+    _, perimeter, area = PlanimeterTest.polygon.TestPoint(1, 1, False, True)
+    self.assertTrue(area == 0)
+    self.assertTrue(perimeter == 0)
+    _, perimeter, area = PlanimeterTest.polygon.TestEdge(90, 1000, False, True)
+    self.assertTrue(math.isnan(area))
+    self.assertTrue(math.isnan(perimeter))
+    PlanimeterTest.polygon.AddPoint(1, 1)
+    _, perimeter, area = PlanimeterTest.polygon.Compute(False, True)
+    self.assertTrue(area == 0)
+    self.assertTrue(perimeter == 0)
+    PlanimeterTest.polyline.Clear()
+    _, perimeter, area = PlanimeterTest.polyline.Compute(False, True)
+    self.assertTrue(perimeter == 0)
+    _, perimeter, area = PlanimeterTest.polyline.TestPoint(1, 1, False, True)
+    self.assertTrue(perimeter == 0)
+    _, perimeter, area = PlanimeterTest.polyline.TestEdge(90, 1000, False, True)
+    self.assertTrue(math.isnan(perimeter))
+    PlanimeterTest.polyline.AddPoint(1, 1)
+    _, perimeter, area = PlanimeterTest.polyline.Compute(False, True)
+    self.assertTrue(perimeter == 0)
+    PlanimeterTest.polygon.AddPoint(1, 1)
+    _, perimeter, area = PlanimeterTest.polyline.TestEdge(90, 1000, False, True)
+    self.assertAlmostEqual(perimeter, 1000, delta = 1e-10)
+    _, perimeter, area = PlanimeterTest.polyline.TestPoint(2, 2, False, True)
+    self.assertAlmostEqual(perimeter, 156876.149, delta = 0.5e-3)
+
+  def test_Planimeter21(self):
+    """Some tests to add code coverage: multiple circlings of pole (includes
+    Planimeter21 - Planimeter28) + invocations via testpoint and testedge."""
+    lat = 45
+    azi = 39.2144607176828184218
+    s = 8420705.40957178156285
+    r = 39433884866571.4277     # Area for one circuit
+    a0 = 510065621724088.5093   # Ellipsoid area
+    PlanimeterTest.polygon.Clear()
+    PlanimeterTest.polygon.AddPoint(lat,  60)
+    PlanimeterTest.polygon.AddPoint(lat, 180)
+    PlanimeterTest.polygon.AddPoint(lat, -60)
+    PlanimeterTest.polygon.AddPoint(lat,  60)
+    PlanimeterTest.polygon.AddPoint(lat, 180)
+    PlanimeterTest.polygon.AddPoint(lat, -60)
+    for i in [3, 4]:
+      PlanimeterTest.polygon.AddPoint(lat,  60)
+      PlanimeterTest.polygon.AddPoint(lat, 180)
+      _, _, area = PlanimeterTest.polygon.TestPoint(lat, -60, False, True)
+      self.assertAlmostEqual(area,  i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestPoint(lat, -60, False, False)
+      self.assertAlmostEqual(area,  i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestPoint(lat, -60, True, True)
+      self.assertAlmostEqual(area, -i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestPoint(lat, -60, True, False)
+      self.assertAlmostEqual(area, -i*r + a0, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestEdge(azi, s, False, True)
+      self.assertAlmostEqual(area,  i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestEdge(azi, s, False, False)
+      self.assertAlmostEqual(area,  i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestEdge(azi, s, True, True)
+      self.assertAlmostEqual(area, -i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.TestEdge(azi, s, True, False)
+      self.assertAlmostEqual(area, -i*r + a0, delta = 0.5)
+      PlanimeterTest.polygon.AddPoint(lat, -60)
+      _, _, area = PlanimeterTest.polygon.Compute(False, True)
+      self.assertAlmostEqual(area,  i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.Compute(False, False)
+      self.assertAlmostEqual(area,  i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.Compute(True, True)
+      self.assertAlmostEqual(area, -i*r, delta = 0.5)
+      _, _, area = PlanimeterTest.polygon.Compute(True, False)
+      self.assertAlmostEqual(area, -i*r + a0, delta = 0.5)
+
+  def test_Planimeter29(self):
+    """Check fix to transitdirect vs transit zero handling inconsistency"""
+    PlanimeterTest.polygon.Clear()
+    PlanimeterTest.polygon.AddPoint(0, 0)
+    PlanimeterTest.polygon.AddEdge( 90, 1000)
+    PlanimeterTest.polygon.AddEdge(  0, 1000)
+    PlanimeterTest.polygon.AddEdge(-90, 1000)
+    _, _, area = PlanimeterTest.polygon.Compute(False, True)
+    # The area should be 1e6.  Prior to the fix it was 1e6 - A/2, where
+    # A = ellipsoid area.
+    self.assertAlmostEqual(area, 1000000.0, delta = 0.01)

streamlit-venv/lib/python3.10/site-packages/geographiclib/test/test_sign.py

@@ -0,0 +1,280 @@
+"""Geodesic tests"""
+
+import unittest
+import math
+import sys
+
+from geographiclib.geomath import Math
+from geographiclib.geodesic import Geodesic
+
+class SignTest(unittest.TestCase):
+  """Sign test suite"""
+
+  @staticmethod
+  def equiv(x, y):
+    """Test for equivalence"""
+
+    return ( (math.isnan(x) and math.isnan(y)) or
+             (x == y and math.copysign(1.0, x) == math.copysign(1.0, y)) )
+
+  def test_AngRound(self):
+    """Test special cases for AngRound"""
+    eps = sys.float_info.epsilon
+    self.assertTrue(SignTest.equiv(Math.AngRound(-eps/32), -eps/32))
+    self.assertTrue(SignTest.equiv(Math.AngRound(-eps/64), -0.0   ))
+    self.assertTrue(SignTest.equiv(Math.AngRound(-  0.0 ), -0.0   ))
+    self.assertTrue(SignTest.equiv(Math.AngRound(   0.0 ), +0.0   ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( eps/64), +0.0   ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( eps/32), +eps/32))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-2*eps)/64), (1-2*eps)/64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps  )/64),  1.0     /64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/2)/64),  1.0     /64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/4)/64),  1.0     /64))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1.0     /64),  1.0     /64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps/2)/64),  1.0     /64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps  )/64),  1.0     /64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+2*eps)/64), (1+2*eps)/64))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps  )/32), (1-eps  )/32))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/2)/32),  1.0     /32))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/4)/32),  1.0     /32))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1.0     /32),  1.0     /32))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps/2)/32),  1.0     /32))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps  )/32), (1+eps  )/32))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps  )/16), (1-eps  )/16))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/2)/16), (1-eps/2)/16))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/4)/16),  1.0     /16))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1.0     /16),  1.0     /16))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps/4)/16),  1.0     /16))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps/2)/16),  1.0     /16))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps  )/16), (1+eps  )/16))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps  )/ 8), (1-eps  )/ 8))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/2)/ 8), (1-eps/2)/ 8))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1-eps/4)/ 8),  1.0     / 8))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps/2)/ 8),  1.0     / 8))
+    self.assertTrue(SignTest.equiv(Math.AngRound((1+eps  )/ 8), (1+eps  )/ 8))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1-eps      ),  1-eps      ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1-eps/2    ),  1-eps/2    ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1-eps/4    ),  1          ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1.0        ),  1          ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1+eps/4    ),  1          ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1+eps/2    ),  1          ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 1+eps      ),  1+  eps    ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 90.0-64*eps),  90-64*eps  ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 90.0-32*eps),  90         ))
+    self.assertTrue(SignTest.equiv(Math.AngRound( 90.0       ),  90         ))
+
+  def test_sincosd(self):
+    """Test special cases for sincosd"""
+    inf = math.inf
+    nan = math.nan
+    s, c = Math.sincosd(-  inf)
+    self.assertTrue(SignTest.equiv(s,  nan) and SignTest.equiv(c,  nan))
+    s, c = Math.sincosd(-810.0)
+    self.assertTrue(SignTest.equiv(s, -1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(-720.0)
+    self.assertTrue(SignTest.equiv(s, -0.0) and SignTest.equiv(c, +1.0))
+    s, c = Math.sincosd(-630.0)
+    self.assertTrue(SignTest.equiv(s, +1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(-540.0)
+    self.assertTrue(SignTest.equiv(s, -0.0) and SignTest.equiv(c, -1.0))
+    s, c = Math.sincosd(-450.0)
+    self.assertTrue(SignTest.equiv(s, -1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(-360.0)
+    self.assertTrue(SignTest.equiv(s, -0.0) and SignTest.equiv(c, +1.0))
+    s, c = Math.sincosd(-270.0)
+    self.assertTrue(SignTest.equiv(s, +1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(-180.0)
+    self.assertTrue(SignTest.equiv(s, -0.0) and SignTest.equiv(c, -1.0))
+    s, c = Math.sincosd(- 90.0)
+    self.assertTrue(SignTest.equiv(s, -1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(-  0.0)
+    self.assertTrue(SignTest.equiv(s, -0.0) and SignTest.equiv(c, +1.0))
+    s, c = Math.sincosd(+  0.0)
+    self.assertTrue(SignTest.equiv(s, +0.0) and SignTest.equiv(c, +1.0))
+    s, c = Math.sincosd(+ 90.0)
+    self.assertTrue(SignTest.equiv(s, +1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(+180.0)
+    self.assertTrue(SignTest.equiv(s, +0.0) and SignTest.equiv(c, -1.0))
+    s, c = Math.sincosd(+270.0)
+    self.assertTrue(SignTest.equiv(s, -1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(+360.0)
+    self.assertTrue(SignTest.equiv(s, +0.0) and SignTest.equiv(c, +1.0))
+    s, c = Math.sincosd(+450.0)
+    self.assertTrue(SignTest.equiv(s, +1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(+540.0)
+    self.assertTrue(SignTest.equiv(s, +0.0) and SignTest.equiv(c, -1.0))
+    s, c = Math.sincosd(+630.0)
+    self.assertTrue(SignTest.equiv(s, -1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(+720.0)
+    self.assertTrue(SignTest.equiv(s, +0.0) and SignTest.equiv(c, +1.0))
+    s, c = Math.sincosd(+810.0)
+    self.assertTrue(SignTest.equiv(s, +1.0) and SignTest.equiv(c, +0.0))
+    s, c = Math.sincosd(+  inf)
+    self.assertTrue(SignTest.equiv(s,  nan) and SignTest.equiv(c,  nan))
+    s, c = Math.sincosd(  nan)
+    self.assertTrue(SignTest.equiv(s,  nan) and SignTest.equiv(c,  nan))
+
+  def test_sincosd2(self):
+    """Test accuracy of sincosd"""
+    s1, c1 = Math.sincosd(         9.0)
+    s2, c2 = Math.sincosd(        81.0)
+    s3, c3 = Math.sincosd(-123456789.0)
+    self.assertTrue(SignTest.equiv(s1, c2))
+    self.assertTrue(SignTest.equiv(s1, s3))
+    self.assertTrue(SignTest.equiv(c1, s2))
+    self.assertTrue(SignTest.equiv(c1,-c3))
+
+  def test_atan2d(self):
+    """Test special cases for atan2d"""
+    inf = math.inf
+    nan = math.nan
+    self.assertTrue(SignTest.equiv(Math.atan2d(+0.0 , -0.0 ), +180))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-0.0 , -0.0 ), -180))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+0.0 , +0.0 ), +0.0))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-0.0 , +0.0 ), -0.0))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+0.0 , -1.0 ), +180))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-0.0 , -1.0 ), -180))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+0.0 , +1.0 ), +0.0))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-0.0 , +1.0 ), -0.0))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-1.0 , +0.0 ),  -90))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-1.0 , -0.0 ),  -90))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+1.0 , +0.0 ),  +90))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+1.0 , -0.0 ),  +90))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+1.0 ,  -inf), +180))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-1.0 ,  -inf), -180))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+1.0 ,  +inf), +0.0))
+    self.assertTrue(SignTest.equiv(Math.atan2d(-1.0 ,  +inf), -0.0))
+    self.assertTrue(SignTest.equiv(Math.atan2d( +inf, +1.0 ),  +90))
+    self.assertTrue(SignTest.equiv(Math.atan2d( +inf, -1.0 ),  +90))
+    self.assertTrue(SignTest.equiv(Math.atan2d( -inf, +1.0 ),  -90))
+    self.assertTrue(SignTest.equiv(Math.atan2d( -inf, -1.0 ),  -90))
+    self.assertTrue(SignTest.equiv(Math.atan2d( +inf,  -inf), +135))
+    self.assertTrue(SignTest.equiv(Math.atan2d( -inf,  -inf), -135))
+    self.assertTrue(SignTest.equiv(Math.atan2d( +inf,  +inf),  +45))
+    self.assertTrue(SignTest.equiv(Math.atan2d( -inf,  +inf),  -45))
+    self.assertTrue(SignTest.equiv(Math.atan2d(  nan, +1.0 ),  nan))
+    self.assertTrue(SignTest.equiv(Math.atan2d(+1.0 ,   nan),  nan))
+
+  def test_atan2d2(self):
+    """Test accuracy of atan2d"""
+    s = 7e-16
+    self.assertEqual(Math.atan2d(s, -1.0), 180 - Math.atan2d(s, 1.0))
+
+  def test_sum(self):
+    """Test special cases of sum"""
+    s,_ = Math.sum(+9.0, -9.0); self.assertTrue(SignTest.equiv(s, +0.0))
+    s,_ = Math.sum(-9.0, +9.0); self.assertTrue(SignTest.equiv(s, +0.0))
+    s,_ = Math.sum(-0.0, +0.0); self.assertTrue(SignTest.equiv(s, +0.0))
+    s,_ = Math.sum(+0.0, -0.0); self.assertTrue(SignTest.equiv(s, +0.0))
+    s,_ = Math.sum(-0.0, -0.0); self.assertTrue(SignTest.equiv(s, -0.0))
+    s,_ = Math.sum(+0.0, +0.0); self.assertTrue(SignTest.equiv(s, +0.0))
+
+  def test_AngNormalize(self):
+    """Test special cases of AngNormalize"""
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(-900.0), -180))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(-720.0), -0.0))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(-540.0), -180))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(-360.0), -0.0))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(-180.0), -180))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(  -0.0), -0.0))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize(  +0.0), +0.0))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize( 180.0), +180))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize( 360.0), +0.0))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize( 540.0), +180))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize( 720.0), +0.0))
+    self.assertTrue(SignTest.equiv(Math.AngNormalize( 900.0), +180))
+
+  def test_AngDiff(self):
+    """Test special cases of AngDiff"""
+    eps = sys.float_info.epsilon
+    s,_ = Math.AngDiff(+  0.0,+  0.0); self.assertTrue(SignTest.equiv(s,+0.0 ))
+    s,_ = Math.AngDiff(+  0.0,-  0.0); self.assertTrue(SignTest.equiv(s,-0.0 ))
+    s,_ = Math.AngDiff(-  0.0,+  0.0); self.assertTrue(SignTest.equiv(s,+0.0 ))
+    s,_ = Math.AngDiff(-  0.0,-  0.0); self.assertTrue(SignTest.equiv(s,+0.0 ))
+    s,_ = Math.AngDiff(+  5.0,+365.0); self.assertTrue(SignTest.equiv(s,+0.0 ))
+    s,_ = Math.AngDiff(+365.0,+  5.0); self.assertTrue(SignTest.equiv(s,-0.0 ))
+    s,_ = Math.AngDiff(+  5.0,+185.0); self.assertTrue(SignTest.equiv(s,+180.0))
+    s,_ = Math.AngDiff(+185.0,+  5.0); self.assertTrue(SignTest.equiv(s,-180.0))
+    s,_ = Math.AngDiff( +eps ,+180.0); self.assertTrue(SignTest.equiv(s,+180.0))
+    s,_ = Math.AngDiff( -eps ,+180.0); self.assertTrue(SignTest.equiv(s,-180.0))
+    s,_ = Math.AngDiff( +eps ,-180.0); self.assertTrue(SignTest.equiv(s,+180.0))
+    s,_ = Math.AngDiff( -eps ,-180.0); self.assertTrue(SignTest.equiv(s,-180.0))
+
+  def test_AngDiff2(self):
+    """Test accuracy of AngDiff"""
+    eps = sys.float_info.epsilon
+    x = 138 + 128 * eps; y = -164; s,_ = Math.AngDiff(x, y)
+    self.assertEqual(s, 58 - 128 * eps)
+
+  def test_equatorial_coincident(self):
+    """
+    azimuth with coincident point on equator
+    """
+    # lat1 lat2 azi1/2
+    C = [
+      [ +0.0, -0.0, 180 ],
+      [ -0.0, +0.0,   0 ]
+    ]
+    for l in C:
+      (lat1, lat2, azi) = l
+      inv = Geodesic.WGS84.Inverse(lat1, 0.0, lat2, 0.0)
+      self.assertTrue(SignTest.equiv(inv["azi1"], azi))
+      self.assertTrue(SignTest.equiv(inv["azi2"], azi))
+
+  def test_equatorial_NS(self):
+    """Does the nearly antipodal equatorial solution go north or south?"""
+    # lat1 lat2 azi1 azi2
+    C = [
+      [ +0.0, +0.0,  56, 124],
+      [ -0.0, -0.0, 124,  56]
+    ]
+    for l in C:
+      (lat1, lat2, azi1, azi2) = l
+      inv = Geodesic.WGS84.Inverse(lat1, 0.0, lat2, 179.5)
+      self.assertAlmostEqual(inv["azi1"], azi1, delta = 1)
+      self.assertAlmostEqual(inv["azi2"], azi2, delta = 1)
+
+  def test_antipodal(self):
+    """How does the exact antipodal equatorial path go N/S + E/W"""
+    # lat1 lat2 lon2 azi1 azi2
+    C = [
+      [ +0.0, +0.0, +180,   +0.0, +180],
+      [ -0.0, -0.0, +180, +180,   +0.0],
+      [ +0.0, +0.0, -180,   -0.0, -180],
+      [ -0.0, -0.0, -180, -180,   -0.0]
+    ]
+    for l in C:
+      (lat1, lat2, lon2, azi1, azi2) = l
+      inv = Geodesic.WGS84.Inverse(lat1, 0.0, lat2, lon2)
+      self.assertTrue(SignTest.equiv(inv["azi1"], azi1))
+      self.assertTrue(SignTest.equiv(inv["azi2"], azi2))
+
+  def test_antipodal_prolate(self):
+    """Antipodal points on the equator with prolate ellipsoid"""
+    # lon2 azi1/2
+    C = [
+      [ +180, +90 ],
+      [ -180, -90 ]
+    ]
+    geod = Geodesic(6.4e6, -1/300.0)
+    for l in C:
+      (lon2, azi) = l
+      inv = geod.Inverse(0.0, 0.0, 0.0, lon2)
+      self.assertTrue(SignTest.equiv(inv["azi1"], azi))
+      self.assertTrue(SignTest.equiv(inv["azi2"], azi))
+
+  def test_azimuth_0_180(self):
+    """azimuths = +/-0 and +/-180 for the direct problem"""
+    # azi1, lon2, azi2
+    C = [
+      [ +0.0, +180, +180 ],
+      [ -0.0, -180, -180 ],
+      [ +180, +180, +0.0 ],
+      [ -180, -180, -0.0 ]
+    ]
+    for l in C:
+      (azi1, lon2, azi2) = l
+      direct = Geodesic.WGS84.Direct(0.0, 0.0, azi1, 15e6,
+                                     Geodesic.STANDARD | Geodesic.LONG_UNROLL)
+      self.assertTrue(SignTest.equiv(direct["lon2"], lon2))
+      self.assertTrue(SignTest.equiv(direct["azi2"], azi2))