From 14f47667f24a0277852bd40fe8efa85ac77b7ba5 Mon Sep 17 00:00:00 2001
From: Julien Deniau <julien.deniau@gmail.com>
Date: Fri, 15 Jul 2022 12:01:55 +0200
Subject: [PATCH] Try to improve inverse_haversine

---
 haversine/haversine.py          | 57 +++++++++++++++++++++++++++++----
 tests/test_inverse_haversine.py | 25 +++++++++++++--
 2 files changed, 72 insertions(+), 10 deletions(-)

diff --git a/haversine/haversine.py b/haversine/haversine.py
index 8140ef9..eabb7a5 100755
--- a/haversine/haversine.py
+++ b/haversine/haversine.py
@@ -199,18 +199,61 @@ def haversine_vector(array1, array2, unit=Unit.KILOMETERS, comb=False, normalize
     return 2 * get_avg_earth_radius(unit) * numpy.arcsin(numpy.sqrt(d))
 
 
-def inverse_haversine(point, distance, direction: Union[Direction, float], unit=Unit.KILOMETERS):
+def _get_cos_direction(direction: Union[Direction, float]) -> float:
+    """
+    Returns the cosine of the direction. Transform special case to avoid mathematical imprecision of computers
+    """
+    brng = direction.value if isinstance(direction, Direction) else direction
+
+    if brng == 0:  # North
+        return 1
+    elif brng == 0.5 * pi:  # East
+        return 0
+    elif brng == pi:  # South
+        return -1
+    elif brng == 1.5 * pi:  # West
+        return 0
+    else:
+        return cos(brng)
+
+
+def _get_sin_direction(direction: Union[Direction, float]) -> float:
+    """
+    Returns the sine of the direction. Transform special case to avoid mathematical imprecision of computers
+    """
+    brng = direction.value if isinstance(direction, Direction) else direction
+
+    if brng == 0:  # North
+        return 0
+    elif brng == 0.5 * pi:  # East
+        return 1
+    elif brng == pi:  # South
+        return 0
+    elif brng == 1.5 * pi:  # West
+        return -1
+    else:
+        return sin(brng)
+
 
+def inverse_haversine(point, distance, direction: Union[Direction, float], unit=Unit.KILOMETERS):
     lat, lng = point
     lat, lng = map(radians, (lat, lng))
-    d = distance
     r = get_avg_earth_radius(unit)
-    brng = direction.value if isinstance(direction, Direction) else direction
 
-    return_lat = asin(sin(lat) * cos(d / r) + cos(lat)
-                      * sin(d / r) * cos(brng))
-    return_lng = lng + atan2(sin(brng) * sin(d / r) *
-                             cos(lat), cos(d / r) - sin(lat) * sin(return_lat))
+    cos_brng = _get_cos_direction(direction)
+    sin_brng = _get_sin_direction(direction)
+    # print(point, distance, r)
+
+    cos_d_r = cos(distance / r)
+    sin_d_r = sin(distance / r)
+    cos_lat = cos(lat)
+    sin_lat = sin(lat)
+
+    # print(sin_lat, cos_d_r, sin_lat * cos_d_r, cos_lat * sin_d_r * cos_brng)
+    return_lat = asin(sin_lat * cos_d_r + cos_lat * sin_d_r * cos_brng)
+    return_lng = lng + atan2(sin_brng * sin_d_r * cos_lat,
+                             cos_d_r - sin_lat * sin(return_lat))
+    # print(degrees(return_lat))
 
     return_lat, return_lng = map(degrees, (return_lat, return_lng))
     return return_lat, return_lng
diff --git a/tests/test_inverse_haversine.py b/tests/test_inverse_haversine.py
index 34cb6fb..c7023eb 100755
--- a/tests/test_inverse_haversine.py
+++ b/tests/test_inverse_haversine.py
@@ -5,6 +5,7 @@
 
 from tests.geo_ressources import LYON, PARIS, NEW_YORK, LONDON
 
+
 @pytest.mark.parametrize(
     "point, dir, dist, result",
     [
@@ -18,7 +19,8 @@
     ],
 )
 def test_inverse_kilometers(point, dir, dist, result):
-    assert isclose(inverse_haversine(point, dist, dir), result, rtol=1e-5).all()
+    assert isclose(inverse_haversine(point, dist, dir),
+                   result, rtol=1e-5).all()
 
 
 @pytest.mark.parametrize(
@@ -36,8 +38,25 @@ def test_back_and_forth(point, direction, distance, unit):
 
 
 def test_inverse_miles():
-    assert isclose(inverse_haversine(PARIS, 50, Direction.NORTH, unit=Unit.MILES), (49.5803579218996, 2.3508), rtol=1e-5).all()
+    assert isclose(inverse_haversine(PARIS, 50, Direction.NORTH,
+                   unit=Unit.MILES), (49.5803579218996, 2.3508), rtol=1e-5).all()
 
 
 def test_nautical_inverse_miles():
-    assert isclose(inverse_haversine(PARIS, 10, Direction.SOUTH, unit=Unit.NAUTICAL_MILES), (48.69014586638915, 2.3508), rtol=1e-5).all()
+    assert isclose(inverse_haversine(PARIS, 10, Direction.SOUTH,
+                   unit=Unit.NAUTICAL_MILES), (48.69014586638915, 2.3508), rtol=1e-5).all()
+
+
+@pytest.mark.parametrize(
+    "point, direction, distance, unit, expected",
+    [
+        (PARIS, Direction.WEST, 3000, Unit.KILOMETERS, (PARIS[0], 0)),
+        # (LYON, Direction.WEST, 3000, Unit.KILOMETERS, (LYON[0], 0)),
+    ],
+)
+def test_explicit_cardinal_points(point, direction, distance, unit, expected):
+    """
+    Test going north/south should not alter latitude and going east/west should not alter longitude
+    """
+    assert inverse_haversine(point, distance, direction, unit)[
+        0] == expected[0]