fix: use geodetic latitudes in haversine distance formula

The implementation was incorrectly using reduced latitudes (via a
flattening factor from WGS84 ellipsoid constants) instead of raw
geodetic latitudes. Reduced latitudes are appropriate for ellipsoidal
models like Lambert's formula, but the Haversine formula operates on
a sphere and should use geodetic latitudes directly.

Changes:
- Use radians(lat) directly instead of computing reduced latitudes
  with atan((1 - flattening) * tan(radians(lat)))
- Replace equatorial radius (6378137m) with mean Earth radius
  (6371000m) for better spherical approximation
- Remove unused WGS84 ellipsoid constants (AXIS_A, AXIS_B)
- Remove unused imports (atan, tan)
- Add edge case and cross-continental doctests

Fixes #11308

Author: VibhorGautam

geodesy/haversine_distance.py

@@ -1,13 +1,11 @@
-from math import asin, atan, cos, radians, sin, sqrt, tan
+from math import asin, cos, radians, sin, sqrt
 
-AXIS_A = 6378137.0
-AXIS_B = 6356752.314245
-RADIUS = 6378137
+EARTH_RADIUS = 6371000
 
 
 def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
     """
-    Calculate great circle distance between two points in a sphere,
+    Calculate great circle distance between two points on a sphere,
     given longitudes and latitudes https://en.wikipedia.org/wiki/Haversine_formula
 
     We know that the globe is "sort of" spherical, so a path between two points
@@ -21,8 +19,10 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl
     computation like Haversine can be handy for shorter range distances.
 
     Args:
-        * `lat1`, `lon1`: latitude and longitude of coordinate 1
-        * `lat2`, `lon2`: latitude and longitude of coordinate 2
+        lat1: latitude of coordinate 1 in degrees
+        lon1: longitude of coordinate 1 in degrees
+        lat2: latitude of coordinate 2 in degrees
+        lon2: longitude of coordinate 2 in degrees
     Returns:
         geographical distance between two points in metres
 
@@ -31,25 +31,39 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl
     >>> SAN_FRANCISCO = point_2d(37.774856, -122.424227)
     >>> YOSEMITE = point_2d(37.864742, -119.537521)
     >>> f"{haversine_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters"
-    '254,352 meters'
+    '253,748 meters'
+    >>> NEW_YORK = point_2d(40.712776, -74.005974)
+    >>> LOS_ANGELES = point_2d(34.052235, -118.243683)
+    >>> f"{haversine_distance(*NEW_YORK, *LOS_ANGELES):0,.0f} meters"
+    '3,935,746 meters'
+    >>> LONDON = point_2d(51.507351, -0.127758)
+    >>> PARIS = point_2d(48.856614, 2.352222)
+    >>> f"{haversine_distance(*LONDON, *PARIS):0,.0f} meters"
+    '343,549 meters'
+    >>> haversine_distance(0, 0, 0, 0)
+    0.0
+    >>> from math import isclose
+    >>> quarter_equator = haversine_distance(0, 0, 0, 90)
+    >>> isclose(quarter_equator, 10_007_543, rel_tol=1e-3)
+    True
     """
-    # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System
-    # Distance in metres(m)
-    # Equation parameters
-    # Equation https://en.wikipedia.org/wiki/Haversine_formula#Formulation
-    flattening = (AXIS_A - AXIS_B) / AXIS_A
-    phi_1 = atan((1 - flattening) * tan(radians(lat1)))
-    phi_2 = atan((1 - flattening) * tan(radians(lat2)))
+    # Convert geodetic coordinates from degrees to radians.
+    # The Haversine formula operates on a sphere, so we use the raw geodetic
+    # latitudes directly rather than reduced latitudes (which apply to
+    # ellipsoidal models like Lambert's formula).
+    # Reference: https://en.wikipedia.org/wiki/Haversine_formula#Formulation
+    phi_1 = radians(lat1)
+    phi_2 = radians(lat2)
     lambda_1 = radians(lon1)
     lambda_2 = radians(lon2)
-    # Equation
+
+    # Haversine equation
     sin_sq_phi = sin((phi_2 - phi_1) / 2)
     sin_sq_lambda = sin((lambda_2 - lambda_1) / 2)
-    # Square both values
     sin_sq_phi *= sin_sq_phi
     sin_sq_lambda *= sin_sq_lambda
     h_value = sqrt(sin_sq_phi + (cos(phi_1) * cos(phi_2) * sin_sq_lambda))
-    return 2 * RADIUS * asin(h_value)
+    return 2 * EARTH_RADIUS * asin(h_value)
 
 
 if __name__ == "__main__":