bugfix: catch Numpy ArithmeticError

fixes #78

src/pymap3d/__init__.py

@@ -29,7 +29,7 @@
 * Fortran: [Maptran3D](https://github.com/geospace-code/maptran3d)
 """
 
-__version__ = "2.10.0"
+__version__ = "3.0.0"
 
 from .aer import aer2ecef, aer2geodetic, ecef2aer, geodetic2aer
 from .ecef import (
@@ -94,6 +94,10 @@
 ]
 
 try:
+    import numpy as np
+
+    np.seterr(all="raise")
+
     from .aer import aer2eci, eci2aer
     from .azelradec import azel2radec, radec2azel
     from .eci import ecef2eci, eci2ecef

src/pymap3d/ecef.py

@@ -12,7 +12,7 @@
 from math import pi
 
 from .ellipsoid import Ellipsoid
-from .mathfun import atan, atan2, cos, degrees, hypot, radians, sin, sqrt, tan
+from .mathfun import atan, atan2, cos, degrees, hypot, isclose, radians, sin, sqrt, tan
 from .utils import sanitize
 
 __all__ = [
@@ -70,8 +70,8 @@ def geodetic2ecef(
         lon = radians(lon)
 
     # radius of curvature of the prime vertical section
-    N = ell.semimajor_axis**2 / sqrt(
-        ell.semimajor_axis**2 * cos(lat) ** 2 + ell.semiminor_axis**2 * sin(lat) ** 2
+    N = ell.semimajor_axis**2 / hypot(
+        ell.semimajor_axis * cos(lat), ell.semiminor_axis * sin(lat)
     )
     # Compute cartesian (geocentric) coordinates given (curvilinear) geodetic coordinates.
     x = (N + alt) * cos(lat) * cos(lon)
@@ -133,7 +133,7 @@ def ecef2geodetic(
     E = sqrt(ell.semimajor_axis**2 - ell.semiminor_axis**2)
 
     # eqn. 4a
-    u = sqrt(0.5 * (r**2 - E**2) + 0.5 * sqrt((r**2 - E**2) ** 2 + 4 * E**2 * z**2))
+    u = sqrt(0.5 * (r**2 - E**2) + 0.5 * hypot(r**2 - E**2, 2 * E * z))
 
     Q = hypot(x, y)
 
@@ -142,8 +142,10 @@ def ecef2geodetic(
     # eqn. 4b
     try:
         Beta = atan(huE / u * z / hypot(x, y))
-    except ZeroDivisionError:
-        if z >= 0:
+    except ArithmeticError:
+        if isclose(z, 0):
+            Beta = 0
+        elif z > 0:
             Beta = pi / 2
         else:
             Beta = -pi / 2

src/pymap3d/latitude.py

@@ -344,7 +344,7 @@ def geodetic2conformal(geodetic_lat, ell: Ellipsoid = None, deg: bool = True):
     #  compute conformal latitudes with correction for points at +90
     try:
         conformal_lat = 2 * atan(sqrt((f4 / f3) * ((f1 / f2) ** e))) - (pi / 2)
-    except ZeroDivisionError:
+    except ArithmeticError:
         conformal_lat = pi / 2
 
     return degrees(conformal_lat) if deg else conformal_lat

src/pymap3d/mathfun.py

@@ -15,6 +15,7 @@
         exp,
         hypot,
         inf,
+        isclose,
         isnan,
         log,
         power,
@@ -36,6 +37,7 @@
         exp,
         hypot,
         inf,
+        isclose,
         isnan,
         log,
         radians,
@@ -48,7 +50,7 @@ def power(x, y):  # type: ignore
         return pow(x, y)
 
     def sign(x) -> float:  # type: ignore
-        """signum function"""
+        """signum"""
         if x < 0:
             y = -1.0
         elif x > 0:
@@ -58,9 +60,12 @@ def sign(x) -> float:  # type: ignore
 
         return y
 
-    def cbrt(x) -> float:  # type: ignore
-        """math.cbrt was added in Python 3.11"""
-        return x ** (1 / 3)
+    try:
+        import math.cbrt as cbrt  # type: ignore
+    except ImportError:
+
+        def cbrt(x) -> float:  # type: ignore
+            return x ** (1 / 3)
 
 
 __all__ = [
@@ -75,6 +80,7 @@ def cbrt(x) -> float:  # type: ignore
     "exp",
     "hypot",
     "inf",
+    "isclose",
     "isnan",
     "log",
     "power",

src/pymap3d/tests/test_geodetic.py

@@ -14,9 +14,21 @@
 B = ELL.semiminor_axis
 
 xyzlla = [
+    ((A, 0, 0), (0, 0, 0)),
     ((A - 1, 0, 0), (0, 0, -1)),
+    ((A + 1, 0, 0), (0, 0, 1)),
+    ((0.1 * A, 0, 0), (0, 0, -0.9 * A)),
+    ((0.001 * A, 0, 0), (0, 0, -0.999 * A)),
+    ((0, A, 0), (0, 90, 0)),
     ((0, A - 1, 0), (0, 90, -1)),
+    ((0, A + 1, 0), (0, 90, 1)),
+    ((0, 0.1 * A, 0), (0, 90, -0.9 * A)),
+    ((0, 0.001 * A, 0), (0, 90, -0.999 * A)),
+    ((0, 0, B), (90, 0, 0)),
+    ((0, 0, B + 1), (90, 0, 1)),
     ((0, 0, B - 1), (90, 0, -1)),
+    ((0, 0, 0.1 * B), (90, 0, -0.9 * B)),
+    ((0, 0, 0.001 * B), (90, 0, -0.999 * B)),
     ((0, 0, B - 1), (89.999999, 0, -1)),
     ((0, 0, B - 1), (89.99999, 0, -1)),
     ((0, 0, -B + 1), (-90, 0, -1)),

src/pymap3d/vincenty.py

@@ -190,10 +190,10 @@ def vdist(
             alpha[isnan(sinAlpha)] = 0
             alpha[(sinAlpha > 1) | (abs(sinAlpha - 1) < 1e-16)] = pi / 2
 
-        except (ZeroDivisionError, TypeError, ValueError):
+        except (ArithmeticError, TypeError, ValueError):
             try:
                 sinAlpha = cos(U1) * cos(U2) * sin(lamb) / sin(sigma)
-            except ZeroDivisionError:
+            except ArithmeticError:
                 sinAlpha = 0.0
 
             if isnan(sinAlpha):