loxodrome_direct: stability near 90,-90,270,-270 azimuth. Fixes #42

Author: scivision

Examples/lox_stability.py

@@ -28,14 +28,15 @@ def matlab_func(lat1: float, lon1: float, rng: float, az: float) -> tuple[float,
 clat, clon, rng = 40.0, -80.0, 10.0  # arbitrary
 
 for i in range(20):
-    azi = 90.0 + 10.0 ** (-i)
-
-    lat, lon = loxodrome_direct(clat, clon, rng, azi)
-
-    assert lat != lat_last
-    assert lon != lon_last
-
-    lat_matlab, lon_matlab = matlab_func(clat, clon, rng, azi)
-    rstr = f"azimuth: {azi} lat,lon: Python: {lat}, {lon}  Matlab: {lat_matlab}, {lon_matlab}"
-    if not (isclose(lat_matlab, lat, rel_tol=0.005) and isclose(lon_matlab, lon, rel_tol=0.001)):
-        logging.error(rstr)
+    for azi in (90 + 10.0 ** (-i), -90 + 10.0 ** (-i), 270 + 10.0 ** (-i), -270 + 10.0 ** (-i)):
+        lat, lon = loxodrome_direct(clat, clon, rng, azi)
+
+        assert lat != lat_last
+        assert lon != lon_last
+
+        lat_matlab, lon_matlab = matlab_func(clat, clon, rng, azi)
+        rstr = f"azimuth: {azi} lat,lon: Python: {lat}, {lon}  Matlab: {lat_matlab}, {lon_matlab}"
+        if not (
+            isclose(lat_matlab, lat, rel_tol=0.005) and isclose(lon_matlab, lon, rel_tol=0.001)
+        ):
+            logging.error(rstr)

src/pymap3d/enu.py

@@ -47,7 +47,7 @@ def enu2aer(
         slant range [meters]
     """
 
-    # 1 millimeter precision for singularity
+    # 1 millimeter precision for singularity stability
 
     try:
         e[abs(e) < 1e-3] = 0.0

src/pymap3d/lox.py

@@ -9,6 +9,7 @@
     from math import radians, degrees, cos, atan2, tan, pi  # type: ignore
 
 from .ellipsoid import Ellipsoid
+from .utils import sign
 from . import rcurve
 from . import rsphere
 from .latitude import (
@@ -146,15 +147,17 @@ def loxodrome_inverse(
 
     # compute azimuth
     az12 = atan2(dlon, disolat)
-    cosaz12 = cos(az12)
+    aux = abs(cos(az12))
 
     # compute distance along loxodromic curve
-    dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / abs(cos(az12))
+    dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / aux
+
+    # straight east or west
     try:
-        if (abs(cosaz12) < 1e-9).any():
-            dist[abs(cosaz12) < 1e-9] = departure(lon2, lon1, lat1, ell, deg=False)
+        if (aux < 1e-9).any():
+            dist[aux < 1e-9] = departure(lon2, lon1, lat1, ell, deg=False)
     except (AttributeError, TypeError):
-        if abs(cosaz12) < 1e-9:  # straight east or west
+        if aux < 1e-9:
             dist = departure(lon2, lon1, lat1, ell, deg=False)
 
     if deg:
@@ -167,14 +170,15 @@ def loxodrome_direct(
     lat1: float | ndarray,
     lon1: float | ndarray,
     rng: float | ndarray,
-    a12: float | float,
+    a12: float,
     ell: Ellipsoid = None,
     deg: bool = True,
-) -> tuple[ndarray, ndarray]:
+) -> tuple[float | ndarray, float | ndarray]:
     """
     Given starting lat, lon with arclength and azimuth, compute final lat, lon
 
     like Matlab reckon('rh', ...)
+    except that "rng" in meters instead of "arclen" degrees of arc
 
     Parameters
     ----------
@@ -226,14 +230,23 @@ def loxodrome_direct(
     newiso = geodetic2isometric(lat2, ell, deg=False)
     iso = geodetic2isometric(lat1, ell, deg=False)
 
+    # stability near singularities
+    i = abs(cos(a12)) < 1e-6
     dlon = tan(a12) * (newiso - iso)
+
+    try:
+        dlon[i] = sign(dlon[i]) * rng[i] / rcurve.transverse(lat1[i], ell=ell, deg=False)  # type: ignore
+    except (AttributeError, TypeError):
+        if i:  # straight east or west
+            dlon = sign(dlon) * rng / rcurve.transverse(lat1, ell=ell, deg=False)
+
     lon2 = lon1 + dlon
 
     if deg:
         lat2, lon2 = degrees(lat2), degrees(lon2)
 
     try:
-        return lat2.squeeze()[()], lon2.squeeze()[()]
+        return lat2.squeeze()[()], lon2.squeeze()[()]  # type: ignore
     except AttributeError:
         return lat2, lon2
 

src/pymap3d/rcurve.py

@@ -53,7 +53,7 @@ def parallel(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
     Returns
     -------
     radius: float
-        radius of ellipsoid
+        radius of ellipsoid (meters)
     """
 
     if deg:
@@ -111,7 +111,7 @@ def transverse(lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True) ->
     Returns
     -------
     radius: float
-        radius of ellipsoid
+        radius of ellipsoid (meters)
     """
 
     lat, ell = sanitize(lat, ell, deg)

src/pymap3d/tests/test_rcurve.py

@@ -41,4 +41,4 @@ def test_numpy_meridian():
 
 def test_numpy_transverse():
     pytest.importorskip("numpy")
-    assert rcurve.transverse([0, 90]) == approx([A, 6399593.626])
+    assert rcurve.transverse([-90, 0, 90]) == approx([6399593.6258, A, 6399593.6258])