diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py
index 39cd250af965..014a65451508 100644
--- a/geodesy/haversine_distance.py
+++ b/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__":
diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py
index 4805674e51ab..c1b1bbdf388b 100644
--- a/geodesy/lamberts_ellipsoidal_distance.py
+++ b/geodesy/lamberts_ellipsoidal_distance.py
@@ -1,6 +1,6 @@
 from math import atan, cos, radians, sin, tan
 
-from .haversine_distance import haversine_distance
+from .haversine_distance import EARTH_RADIUS, haversine_distance
 
 AXIS_A = 6378137.0
 AXIS_B = 6356752.314245
@@ -39,11 +39,11 @@ def lamberts_ellipsoidal_distance(
     >>> NEW_YORK = point_2d(40.713019, -74.012647)
     >>> VENICE = point_2d(45.443012, 12.313071)
     >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters"
-    '254,351 meters'
+    '254,032 meters'
     >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *NEW_YORK):0,.0f} meters"
-    '4,138,992 meters'
+    '4,133,295 meters'
     >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *VENICE):0,.0f} meters"
-    '9,737,326 meters'
+    '9,719,525 meters'
     """
 
     # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System
@@ -58,7 +58,7 @@ def lamberts_ellipsoidal_distance(
 
     # Compute central angle between two points
     # using haversine theta. sigma =  haversine_distance / equatorial radius
-    sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS
+    sigma = haversine_distance(lat1, lon1, lat2, lon2) / EARTH_RADIUS
 
     # Intermediate P and Q values
     p_value = (b_lat1 + b_lat2) / 2