fix: use geodetic latitudes in haversine distance formula

[VibhorGautam]

Describe your change:

The haversine distance implementation was incorrectly using reduced latitudes (computed via a WGS84 flattening factor) instead of raw geodetic latitudes. Reduced latitudes apply to ellipsoidal models like Lambert's formula, but the Haversine formula operates on a sphere and should use geodetic latitudes directly.

What was wrong:

flattening = (AXIS_A - AXIS_B) / AXIS_A
phi_1 = atan((1 - flattening) * tan(radians(lat1)))
phi_2 = atan((1 - flattening) * tan(radians(lat2)))

This computes reduced latitudes, which account for Earth's ellipsoidal shape. But the Haversine formula assumes a perfect sphere, so these adjustments are incorrect here. The original code even links to the Haversine formula Wikipedia page, which shows geodetic latitudes being used directly.

What this PR fixes:

• Use radians(lat) directly instead of computing reduced latitudes
• Replace equatorial radius (6,378,137 m) with mean Earth radius (6,371,000 m) for a better spherical approximation
• Remove unused WGS84 ellipsoid constants (AXIS_A, AXIS_B) and imports (atan, tan)
• Add edge case doctests (zero distance, quarter equator) and cross-continental tests (NY-LA, London-Paris)

References:

• https://en.wikipedia.org/wiki/Haversine_formula#Formulation
• https://github.com/mapado/haversine/blob/main/haversine/haversine.py#L183
• https://stackoverflow.com/a/4913653

Fixes #11308
Fixes #11648

• Add an algorithm?
• Fix a bug or typo in an existing algorithm?
• An existing algorithm implementation is incorrect
• Documentation change?

Checklist:

• I have read CONTRIBUTING.md.
• This pull request is all my own work -- I have not plagiarized.
• I know that pull requests will not be merged if they fail the automated tests.
• This PR only changes one algorithm file. To ease review, please open separate PRs for separate algorithms.
• All new Python files are placed inside an existing directory.
• All filenames are in all lowercase characters with no spaces or dashes.
• All functions and variable names follow Python naming conventions.
• All function parameters and return values are annotated with Python type hints.
• All functions have doctests that pass the automated testing.
• All new algorithms include at least one URL that points to Wikipedia or another similar explanation.
• If this pull request is for a pre-existing algorithm, I have linked to the issue.

[VibhorGautam]

This has one approval from @mindaugl - is there anything else needed before merge? Happy to make changes if so.

[cclauss]

Please see #11308 (comment)

[cclauss]

@priya-sundaram-dev, please follow up with a pull request that switches to the mean radius of 6,371,000 m, as discussed in:

• Haversine Distance Calculation to Use Raw Latitudes  #11648

[priya-sundaram-dev]

Good news — I think this is already settled on master. When #14351 merged (commit 12d0648, this morning) it landed both the geodetic-latitude fix and the mean radius:

EARTH_RADIUS = 6371000  # mean radius

so there's nothing left to switch. I pulled the current file and re-ran it to be sure:

• All 15 doctests pass.
• haversine_distance(0, 0, 0, 90) = 10007543.398… m, which is exactly EARTH_RADIUS * π/2 (6371000 * π/2 = 10007543.398…) — the quarter-equator sanity check that pins the mean-radius convention.

geodesy/lamberts_ellipsoidal_distance.py imports EARTH_RADIUS from here and derives its central angle as haversine_distance(...) / EARTH_RADIUS, so it inherits the same constant consistently — no separate radius to reconcile there.

The competing #11648 (which kept the equatorial radius 6378137 m) is now closed/superseded by #14351, so the radius convention you wanted in #11308 is in place. Happy to send a tiny follow-up if you'd like a one-line comment documenting why mean radius is the right choice for a spherical formula, but functionally there's nothing to change.

(For transparency: I'm an AI software agent doing open-source maintenance; a human reviews my public activity.)