add example for loxodrome_direct / reckon stability re: #42

Author: scivision

Examples/lox_stability.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+from __future__ import annotations
+import logging
+from pathlib import Path
+from math import isclose, nan
+import numpy as np
+
+from pymap3d.lox import loxodrome_direct
+
+import matlab.engine
+
+cwd = Path(__file__).parent
+eng = None  # don't start Matlab engine over and over when script is interactive
+
+if eng is None:
+    eng = matlab.engine.start_matlab("-nojvm")
+    eng.addpath(eng.genpath(str(cwd)), nargout=0)
+
+if not eng.has_map_toolbox():
+    raise EnvironmentError("Matlab does not have Mapping Toolbox")
+
+
+def matlab_func(lat1, lon1, rng, az) -> tuple[float, float]:
+    """Using Matlab Engine to do same thing as Pymap3d"""
+    return eng.reckon("rh", lat1, lon1, rng, az, eng.wgs84Ellipsoid(), nargout=2)  # type: ignore
+
+
+lat_last, lon_last = nan, nan
+clat, clon, rng = np.array(40.0), np.array(-80.0), np.array(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)

Examples/vdist_stability.py

@@ -18,10 +18,9 @@
     raise EnvironmentError("Matlab does not have Mapping Toolbox")
 
 
-def matlab_func(lat1: float, lon1: float, lat2: float, lon2: float) -> tuple[float, float]:
+def matlab_func(lat1, lon1, lat2, lon2) -> tuple[float, float]:
     """Using Matlab Engine to do same thing as Pymap3d"""
-    ell = eng.wgs84Ellipsoid()  # type: ignore
-    return eng.distance(lat1, lon1, lat2, lon2, ell, nargout=2)  # type: ignore
+    return eng.distance(lat1, lon1, lat2, lon2, eng.wgs84Ellipsoid(), nargout=2)  # type: ignore
 
 
 dlast, alast = nan, nan

src/pymap3d/tests/test_rhumb.py

@@ -77,7 +77,7 @@ def test_numpy_loxodrome_inverse():
 @pytest.mark.parametrize(
     "lat0,lon0,rng,az,lat1,lon1",
     [
-        (40, -80, 10000, 30, 40.077995, -79.9414144),
+        (40, -80, 10000, 30, 40.0000779959676, -79.9999414477481),
         (0, 0, 0, 0, 0, 0),
         (0, 0, 10018754.17, 90, 0, 90),
         (0, 0, 10018754.17, -90, 0, -90),
@@ -86,9 +86,12 @@ def test_numpy_loxodrome_inverse():
     ],
 )
 def test_loxodrome_direct(lat0, lon0, rng, az, lat1, lon1):
+    """
+    reckon('rh', 40, -80, 10, 30, wgs84Ellipsoid)
+    """
     lat2, lon2 = lox.loxodrome_direct(lat0, lon0, rng, az)
-    assert lat2 == approx(lat1, abs=1e-6)
-    assert lon2 == approx(lon1)
+    assert lat2 == approx(lat1, rel=0.005, abs=1e-6)
+    assert lon2 == approx(lon1, rel=0.001)
     assert isinstance(lat2, float)
     assert isinstance(lon2, float)