loxodrome_direct: further correct east-west asymptote

fixes #43

Author: scivision

Examples/lox_stability.py

@@ -2,7 +2,7 @@
 from __future__ import annotations
 import logging
 from pathlib import Path
-from math import isclose, nan
+from math import isclose
 
 from pymap3d.lox import loxodrome_direct
 
@@ -24,16 +24,12 @@ def matlab_func(lat1: float, lon1: float, rng: float, az: float) -> tuple[float,
     return eng.reckon("rh", lat1, lon1, rng, az, eng.wgs84Ellipsoid(), nargout=2)  # type: ignore
 
 
-lat_last, lon_last = nan, nan
-clat, clon, rng = 40.0, -80.0, 10.0  # arbitrary
+clat, clon, rng = 35.0, 140.0, 50000.0  # arbitrary
 
 for i in range(20):
     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 (

src/pymap3d/lox.py

@@ -4,9 +4,11 @@
 import typing
 
 try:
-    from numpy import radians, degrees, cos, arctan2 as atan2, tan, pi, array, broadcast_arrays
+    from numpy import radians, degrees, cos, arctan2 as atan2, tan, array, broadcast_arrays
 except ImportError:
-    from math import radians, degrees, cos, atan2, tan, pi  # type: ignore
+    from math import radians, degrees, cos, atan2, tan  # type: ignore
+
+from math import pi, tau
 
 from .ellipsoid import Ellipsoid
 from .utils import sign
@@ -34,6 +36,9 @@
 ]
 
 
+COS_EPS = 1e-9
+
+
 def meridian_dist(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
     """
     Computes the ground distance on an ellipsoid from the equator to the input latitude.
@@ -159,7 +164,7 @@ def loxodrome_inverse(
     dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / aux
 
     # straight east or west
-    i = aux < 1e-9
+    i = aux < COS_EPS
     try:
         dist[i] = departure(lon2[i], lon1[i], lat1[i], ell, deg=False)
     except (AttributeError, TypeError):
@@ -215,8 +220,10 @@ def loxodrome_direct(
     if deg:
         lat1, lon1, a12 = radians(lat1), radians(lon1), radians(a12)
 
+    a12 = a12 % tau
+
     try:
-        lat1, rng = broadcast_arrays(lat1, rng)
+        lat1, rng, a12 = broadcast_arrays(lat1, rng, a12)
         if (abs(lat1) > pi / 2).any():  # type: ignore
             raise ValueError("-90 <= latitude <= 90")
         if (rng < 0).any():  # type: ignore
@@ -239,14 +246,14 @@ def loxodrome_direct(
     iso = geodetic2isometric(lat1, ell, deg=False)
 
     # stability near singularities
-    i = abs(cos(a12)) < 1e-9
+    i = abs(cos(a12)) < COS_EPS
     dlon = tan(a12) * (newiso - iso)
 
     try:
-        dlon[i] = sign(dlon[i]) * rng[i] / rcurve.transverse(lat1[i], ell=ell, deg=False)  # type: ignore
+        dlon[i] = sign(pi - a12[i]) * rng[i] / rcurve.parallel(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)
+            dlon = sign(pi - a12) * rng / rcurve.parallel(lat1, ell=ell, deg=False)
 
     lon2 = lon1 + dlon
 

src/pymap3d/rcurve.py

@@ -4,9 +4,9 @@
 import typing
 
 try:
-    from numpy import radians, sin, cos, sqrt
+    from numpy import sin, cos, sqrt
 except ImportError:
-    from math import radians, sin, cos, sqrt  # type: ignore
+    from math import sin, cos, sqrt  # type: ignore
 
 from .ellipsoid import Ellipsoid
 from .utils import sanitize
@@ -37,10 +37,12 @@ def geocentric_radius(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool =
     )
 
 
-def parallel(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
+def parallel(lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
     """
     computes the radius of the small circle encompassing the globe at the specified latitude
 
+    like Matlab rcurve('parallel', ...)
+
     Parameters
     ----------
     lat : float
@@ -56,8 +58,7 @@ def parallel(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
         radius of ellipsoid (meters)
     """
 
-    if deg:
-        lat = radians(lat)
+    lat, ell = sanitize(lat, ell, deg)
 
     return cos(lat) * transverse(lat, ell, deg=False)
 
@@ -82,10 +83,7 @@ def meridian(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
         radius of ellipsoid
     """
 
-    if ell is None:
-        ell = Ellipsoid()
-    if deg:
-        lat = radians(lat)
+    lat, ell = sanitize(lat, ell, deg)
 
     f1 = ell.semimajor_axis * (1 - ell.eccentricity ** 2)
     f2 = 1 - (ell.eccentricity * sin(lat)) ** 2

src/pymap3d/tests/test_rhumb.py

@@ -95,6 +95,10 @@ def test_numpy_2d_loxodrome_inverse():
     "lat0,lon0,rng,az,lat1,lon1",
     [
         (40, -80, 10000, 30, 40.0000779959676, -79.9999414477481),
+        (35, 140, 50000, 90, 35, 140.548934481815),
+        (35, 140, 50000, -270, 35, 140.548934481815),
+        (35, 140, 50000, 270, 35, 139.451065518185),
+        (35, 140, 50000, -90, 35, 139.451065518185),
         (0, 0, 0, 0, 0, 0),
         (0, 0, 10018754.17, 90, 0, 90),
         (0, 0, 10018754.17, -90, 0, -90),