loxodrome_inverse: correct singularity handling

Author: scivision

src/pymap3d/latitude.py

@@ -209,12 +209,12 @@ def geodetic2isometric(
     # isometric_lat = log(tan(a2)) + e/2 * log((1-e*sin(geodetic_lat)) / (1+e*sin(geodetic_lat)))
 
     try:
-        isometric_lat[abs(geodetic_lat - pi / 2) <= 1e-9] = inf
-        isometric_lat[abs(-geodetic_lat - pi / 2) <= 1e-9] = -inf
+        isometric_lat[abs(geodetic_lat - pi / 2) <= 1e-9] = inf  # type: ignore
+        isometric_lat[abs(-geodetic_lat - pi / 2) <= 1e-9] = -inf   # type: ignore
     except TypeError:
-        if abs(geodetic_lat - pi / 2) <= 1e-9:
+        if abs(geodetic_lat - pi / 2) <= 1e-9:  # type: ignore
             isometric_lat = inf
-        elif abs(-geodetic_lat - pi / 2) <= 1e-9:
+        elif abs(-geodetic_lat - pi / 2) <= 1e-9:  # type: ignore
             isometric_lat = -inf
 
     return degrees(isometric_lat) if deg else isometric_lat

src/pymap3d/lox.py

@@ -97,6 +97,8 @@ def loxodrome_inverse(
 
     like Matlab distance('rh',...) and azimuth('rh',...)
 
+    If any of inputs lat1,lon1,lat2,lon2 are arrays, all must be arrays of same shape
+
     Parameters
     ----------
 
@@ -153,11 +155,11 @@ def loxodrome_inverse(
     dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / aux
 
     # straight east or west
+    i = aux < 1e-9
     try:
-        if (aux < 1e-9).any():
-            dist[aux < 1e-9] = departure(lon2, lon1, lat1, ell, deg=False)
+        dist[i] = departure(lon2[i], lon1[i], lat1[i], ell, deg=False)
     except (AttributeError, TypeError):
-        if aux < 1e-9:
+        if i:
             dist = departure(lon2, lon1, lat1, ell, deg=False)
 
     if deg:
@@ -180,6 +182,8 @@ def loxodrome_direct(
     like Matlab reckon('rh', ...)
     except that "rng" in meters instead of "arclen" degrees of arc
 
+    If any of inputs lat,lon1,rng are arrays, all must be arrays of same shape
+
     Parameters
     ----------
     lat1 : float
@@ -231,7 +235,7 @@ def loxodrome_direct(
     iso = geodetic2isometric(lat1, ell, deg=False)
 
     # stability near singularities
-    i = abs(cos(a12)) < 1e-6
+    i = abs(cos(a12)) < 1e-9
     dlon = tan(a12) * (newiso - iso)
 
     try:

src/pymap3d/utils.py

@@ -71,10 +71,10 @@ def sanitize(
         lat = radians(lat)
 
     if asarray is not None:
-        if (abs(lat) > pi / 2).any():
+        if (abs(lat) > pi / 2).any():  # type: ignore
             raise ValueError("-pi/2 <= latitude <= pi/2")
     else:
-        if abs(lat) > pi / 2:
+        if abs(lat) > pi / 2:  # type: ignore
             raise ValueError("-pi/2 <= latitude <= pi/2")
 
     return lat, ell