vdist, vreckon: add deg= option (default true)

track2: correct deg=False radians

pure Python linspace()

src/pymap3d/mathfun.py

@@ -2,6 +2,8 @@
 import from Numpy, and if not available fallback to math stdlib
 """
 
+from __future__ import annotations
+
 try:
     from numpy import arcsin as asin
     from numpy import arcsinh as asinh
@@ -17,6 +19,7 @@
         inf,
         isclose,
         isnan,
+        linspace,
         log,
         power,
         radians,
@@ -46,6 +49,14 @@
         tan,
     )
 
+    def linspace(start: float, stop: float, num: int) -> list[float]:  # type: ignore
+        """
+        create a list of "num" evenly spaced numbers using range and increment,
+        including endpoint "stop"
+        """
+        step = (stop - start) / (num - 1)
+        return [start + i * step for i in range(num)]
+
     def power(x, y):  # type: ignore
         return pow(x, y)
 

src/pymap3d/tests/test_matlab_track2.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+"""
+Compare with Matlab Mapping toolbox reckon()
+"""
+
+from __future__ import annotations
+
+import pytest
+from pytest import approx
+
+try:
+    import numpy as np
+    from .matlab_engine import matlab_engine, has_mapping
+except ImportError:
+    pytest.skip("Matlab Engine not found", allow_module_level=True)
+
+
+import pymap3d.vincenty
+
+
+def track2(eng, lat1: float, lon1: float, lat2: float, lon2: float, npts: int, deg: bool) -> tuple:
+    """Using Matlab Engine to do same thing as Pymap3d"""
+    d = "degrees" if deg else "radians"
+
+    lats, lons = eng.track2(
+        "gc", lat1, lon1, lat2, lon2, eng.wgs84Ellipsoid(), d, float(npts), nargout=2
+    )
+    return np.array(lats).squeeze(), np.array(lons).squeeze()
+
+
+@pytest.mark.parametrize("deg", [True, False])
+def test_track2_compare(deg):
+    lat1, lon1 = 0.0, 80.0
+    lat2, lon2 = 0.0, 81.0
+    if not deg:
+        lat1, lon1, lat2, lon2 = np.radians((lat1, lon1, lat2, lon2))
+
+    eng = matlab_engine()
+
+    if not has_mapping(eng):
+        pytest.skip("Matlab Toolbox not found")
+
+    lats, lons = pymap3d.vincenty.track2(lat1, lon1, lat2, lon2, npts=4, deg=deg)
+
+    lats_m, lons_m = track2(eng, lat1, lon1, lat2, lon2, npts=4, deg=deg)
+
+    assert lats == approx(lats_m)
+    assert lons == approx(lons_m)

src/pymap3d/tests/test_vincenty.py

@@ -1,8 +1,25 @@
 import pymap3d.vincenty as vincenty
+
+from math import radians
+
+import pytest
 from pytest import approx
 
 
-def test_track2():
-    lats, lons = vincenty.track2(40, 80, 65, -148, npts=3)
-    assert lats == approx([40, 69.633139886, 65])
-    assert lons == approx([80, 113.06849104, -148])
+@pytest.mark.parametrize("deg", [True, False])
+def test_track2_unit(deg):
+    pytest.importorskip("numpy")
+
+    lat1, lon1 = 0.0, 80.0
+    lat2, lon2 = 0.0, 81.0
+    lat0 = [0.0, 0.0, 0.0, 0.0]
+    lon0 = [80.0, 80.33333, 80.66666, 81.0]
+    if not deg:
+        lat1, lon1, lat2, lon2 = map(radians, (lat1, lon1, lat2, lon2))
+        lat0 = map(radians, lat0)
+        lon0 = map(radians, lon0)
+
+    lats, lons = vincenty.track2(lat1, lon1, lat2, lon2, npts=4, deg=deg)
+
+    assert lats == approx(lat0)
+    assert lons == approx(lon0)

src/pymap3d/tests/test_vincenty_vreckon.py

@@ -1,4 +1,7 @@
+from math import radians
+
 import pymap3d.vincenty as vincenty
+
 import pytest
 from pytest import approx
 
@@ -14,6 +17,7 @@
 az3 = (218.00292856, 225.0011203)
 
 
+@pytest.mark.parametrize("deg", [True, False])
 @pytest.mark.parametrize(
     "lat,lon,srange,az,lato,lono",
     [
@@ -26,8 +30,11 @@
         (0, 0, 2.00375e7, -90, 0, 180),
     ],
 )
-def test_unit(lat, lon, srange, az, lato, lono):
-    lat1, lon1 = vincenty.vreckon(lat, lon, srange, az)
+def test_vreckon_unit(deg, lat, lon, srange, az, lato, lono):
+    if not deg:
+        lat, lon, az, lato, lono = map(radians, (lat, lon, az, lato, lono))
+
+    lat1, lon1 = vincenty.vreckon(lat, lon, srange, az, deg=deg)
 
     assert lat1 == approx(lato)
     assert isinstance(lat1, float)