If you have two GPS coordinates and want the distance between them, it is tempting to treat latitude and longitude like ordinary x/y coordinates and use the Euclidean distance formula.
That shortcut breaks down once the Earth’s curvature matters. Latitude and longitude are positions on the surface of an approximately spherical planet, so the shortest surface path between two points is not a straight line on a flat plane. It is an arc along the sphere.
The haversine formula is a practical way to compute that arc length. Given two lat/lon pairs, it tells you the great-circle distance between them: the shortest distance along the surface of a sphere.
What the haversine formula measures
Imagine two points on the Earth:
- Point 1: latitude \(\phi_1\), longitude \(\lambda_1\)
- Point 2: latitude \(\phi_2\), longitude \(\lambda_2\)
The shortest path constrained to the surface of a sphere lies on a great circle, which is any circle whose center matches the sphere’s center. The equator is a great circle. So is the meridian plane passing through London and New York.
Try dragging the points in the demo below. It shows how the surface arc and the central angle change together.
If we can find the central angle \(\theta\) between the two points, measured from the center of the Earth, then the surface distance is just arc length:
\[d = r\theta\]where \(r\) is the radius of the sphere and \(\theta\) is measured in radians.
This is the standard arc length formula. A full circle has circumference \(2\pi r\) and angle \(2\pi\) radians, so an angle of \(\theta\) covers the fraction \(\theta / 2\pi\) of the circle. The corresponding fraction of the circumference is:
\[d = \frac{\theta}{2\pi} \cdot 2\pi r = r\theta\]So the problem reduces to: given two lat/lon pairs, how do we compute \(\theta\)?
From spherical geometry to haversine
The spherical law of cosines gives the central angle between two points on a sphere:
\[\cos(\theta) = \sin(\phi_1)\sin(\phi_2) + \cos(\phi_1)\cos(\phi_2)\cos(\lambda_2 - \lambda_1)\]This is correct, but in practice it can lose precision for very small distances, where \(\theta\) is close to zero and \(\cos(\theta)\) is close to 1.
The haversine form rewrites the same relationship in a numerically friendlier way.
The haversine function is defined as:
\[\operatorname{hav}(x) = \sin^2\left(\frac{x}{2}\right)\]To get from the spherical law of cosines to the haversine form, we use the cosine difference identity and the half-angle relationships below:
\[\cos(a - b) = \cos(a)\cos(b) + \sin(a)\sin(b)\] \[\cos(x) = 1 - 2\sin^2\left(\frac{x}{2}\right)\] \[\operatorname{hav}(x) = \sin^2\left(\frac{x}{2}\right) = \frac{1 - \cos(x)}{2}\]From the last identity, we can rearrange to get:
\[1 - \cos(x) = 2\operatorname{hav}(x)\]Let:
\[\Delta \phi = \phi_2 - \phi_1\] \[\Delta \lambda = \lambda_2 - \lambda_1\]First apply the cosine difference identity to \(\Delta \phi\):
\[\cos(\Delta \phi) = \cos(\phi_1)\cos(\phi_2) + \sin(\phi_1)\sin(\phi_2)\]Rewrite the spherical law of cosines using \(\Delta \lambda\):
\[\cos(\theta) = \sin(\phi_1)\sin(\phi_2) + \cos(\phi_1)\cos(\phi_2)\cos(\Delta \lambda)\]Now solve the previous identity for \(\sin(\phi_1)\sin(\phi_2)\):
\[\sin(\phi_1)\sin(\phi_2) = \cos(\Delta \phi) - \cos(\phi_1)\cos(\phi_2)\]Substitute that into the spherical law of cosines:
\[\cos(\theta) = \cos(\Delta \phi) - \cos(\phi_1)\cos(\phi_2) + \cos(\phi_1)\cos(\phi_2)\cos(\Delta \lambda)\]Now factor out \(\cos(\phi_1)\cos(\phi_2)\):
\[\cos(\theta) = \cos(\Delta \phi) - \cos(\phi_1)\cos(\phi_2)\left(1 - \cos(\Delta \lambda)\right)\]Next convert everything into haversines by using \(1 - \cos(x) = 2\operatorname{hav}(x)\):
\[1 - \cos(\theta) = 1 - \cos(\Delta \phi) + \cos(\phi_1)\cos(\phi_2)\left(1 - \cos(\Delta \lambda)\right)\] \[2\operatorname{hav}(\theta) = 2\operatorname{hav}(\Delta \phi) + 2\cos(\phi_1)\cos(\phi_2)\operatorname{hav}(\Delta \lambda)\]Divide both sides by 2:
\[\operatorname{hav}(\theta) = \operatorname{hav}(\phi_2 - \phi_1) + \cos(\phi_1)\cos(\phi_2)\operatorname{hav}(\lambda_2 - \lambda_1)\]Now replace each haversine term with \(\operatorname{hav}(x) = \sin^2\left(\frac{x}{2}\right)\):
\[\operatorname{hav}(\theta) = \sin^2\left(\frac{\theta}{2}\right)\] \[\operatorname{hav}(\Delta \phi) = \sin^2\left(\frac{\Delta \phi}{2}\right)\] \[\operatorname{hav}(\Delta \lambda) = \sin^2\left(\frac{\Delta \lambda}{2}\right)\]So:
\[\sin^2\left(\frac{\theta}{2}\right) = \sin^2\left(\frac{\Delta \phi}{2}\right) + \cos(\phi_1)\cos(\phi_2)\sin^2\left(\frac{\Delta \lambda}{2}\right)\]In code, we usually name the right-hand side \(a\):
\[a = \sin^2\left(\frac{\Delta \phi}{2}\right) + \cos(\phi_1)\cos(\phi_2)\sin^2\left(\frac{\Delta \lambda}{2}\right)\]Then:
\[\sin^2\left(\frac{\theta}{2}\right) = a\]Since the central angle on a sphere satisfies \(0 \le \theta \le \pi\), we have \(0 \le \theta/2 \le \pi/2\), so the sine is nonnegative. Taking square roots gives:
\[\sin\left(\frac{\theta}{2}\right) = \sqrt{a}\]Apply inverse sine to recover the central angle:
\[\theta = 2\arcsin(\sqrt{a})\]Many implementations use the equivalent and numerically stable form:
\[c = 2\operatorname{atan2}(\sqrt{a}, \sqrt{1-a})\] \[d = rc\]where \(c = \theta\) is the central angle in radians.
Why half-angles appear
The half-angle terms are not arbitrary. They come from the identity:
\[\operatorname{hav}(x) = \frac{1 - \cos(x)}{2}\]This matters because for small angles, directly subtracting from 1 can amplify floating-point error. Writing the expression in terms of \(\sin^2\left(\frac{x}{2}\right)\) behaves better numerically.
That is one reason the haversine formula became the standard practical choice for many navigation and mapping problems.
Coordinate units matter
Latitude and longitude are usually stored in degrees, but JavaScript, C, Python, and most math libraries expect trigonometric inputs in radians.
So before applying the formula:
\[\text{radians} = \text{degrees} \cdot \frac{\pi}{180}\]If you forget this conversion, the result will be wrong by a large factor.
A practical JavaScript implementation
Here is a compact implementation:
const EARTH_RADIUS_METERS = 6371008.8;
function toRadians(degrees) {
return degrees * Math.PI / 180;
}
function haversineDistance(
lat1,
lon1,
lat2,
lon2,
radius = EARTH_RADIUS_METERS
) {
const phi1 = toRadians(lat1);
const phi2 = toRadians(lat2);
const dPhi = toRadians(lat2 - lat1);
const dLambda = toRadians(lon2 - lon1);
const sinHalfDPhi = Math.sin(dPhi / 2);
const sinHalfDLambda = Math.sin(dLambda / 2);
const a =
sinHalfDPhi * sinHalfDPhi +
Math.cos(phi1) * Math.cos(phi2) *
sinHalfDLambda * sinHalfDLambda;
const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return radius * c;
}
This returns the distance in meters when using the default Earth radius \(r = 6371008.8\) meters.
Example:
const london = { lat: 51.5074, lon: -0.1278 };
const newYork = { lat: 40.7128, lon: -74.0060 };
const distanceMeters = haversineDistance(
london.lat,
london.lon,
newYork.lat,
newYork.lon
);
console.log(distanceMeters / 1000); // about 5570 km
Interpreting the terms
The formula has a clean geometric interpretation:
- \(\sin^2\left(\frac{\Delta \phi}{2}\right)\) captures north-south separation
- \(\sin^2\left(\frac{\Delta \lambda}{2}\right)\) captures east-west separation
- \(\cos(\phi_1)\cos(\phi_2)\) scales longitude separation by latitude
That last term is important. Lines of longitude converge toward the poles, so one degree of longitude does not represent a fixed physical distance everywhere on Earth. Near the equator it spans a large distance. Near the poles it shrinks dramatically.
Why not use plain Euclidean distance?
Suppose two points differ by 1 degree of longitude.
At the equator, that corresponds to roughly 111 km. Near latitude of 60 degrees, it is roughly half of that. Near the poles, it approaches zero. A flat 2D distance formula ignores this geometry entirely, so errors become larger as distances grow or as you move away from the equator.
Projecting to a local plane
For small regions, a common approach is to approximate the Earth’s surface as flat and work in a local Cartesian coordinate system. The idea is to convert latitude and longitude offsets into metric distances at the reference latitude of the region.
Given a reference point \((\phi_0, \lambda_0)\) and a nearby point \((\phi, \lambda)\), the offsets in meters are approximately:
\[\Delta x = (\lambda - \lambda_0) \cdot \cos(\phi_0) \cdot M\] \[\Delta y = (\phi - \phi_0) \cdot M\]where \(M \approx 111{,}320\) meters per degree (one degree of arc along a great circle). The \(\cos(\phi_0)\) factor accounts for the fact that longitude lines converge toward the poles: the east-west distance per degree of longitude shrinks as latitude increases.
The Euclidean distance in the local plane is then:
\[d \approx \sqrt{\Delta x^2 + \Delta y^2}\]This approximation is called the equirectangular projection (or plate carree projection). It is fast and simple, but the error grows with distance from the reference point and with latitude. The flat-earth assumption breaks down once the angular separation between the points is large enough that curvature matters.
As a rough guideline:
| Separation | Typical equirectangular error |
|---|---|
| < 10 km | < 0.01 % |
| ~ 100 km | ~ 0.3 % |
| ~ 1000 km | ~ 3 % or more |
Beyond a few tens of kilometers, or near the poles where the \(\cos(\phi_0)\) factor changes rapidly over the region, haversine is the safer default.
Accuracy and limitations
The haversine formula assumes the Earth is a perfect sphere. It is not. The real Earth is better modeled as an oblate spheroid, slightly flattened at the poles. The polar radius is about 6357 km and the equatorial radius is about 6378 km, a difference of roughly 21 km (0.3 %).
Haversine uses a mean spherical radius of approximately 6371 km. The mismatch between sphere and spheroid introduces a systematic error that varies with the direction of travel.
The worst-case error of haversine against the WGS84 ellipsoid is roughly 0.5 %, which corresponds to about 5 meters per kilometer of distance. The error is largest for long paths that run diagonally between equatorial and polar latitudes; it is near zero for paths along the equator or along a meridian, because those paths happen to lie on circles of the spheroid.
As a practical table:
| Distance | Max haversine error (approx.) |
|---|---|
| 1 km | ~ 5 m |
| 10 km | ~ 50 m |
| 100 km | ~ 500 m |
| 1000 km | ~ 5 km |
These are worst-case estimates. Typical real-world paths have errors closer to half of these figures.
For comparison, the equirectangular approximation over the same 1000 km would accumulate errors of 30 km or more. Haversine is therefore a large improvement over flat-Earth approximations, and accurate enough for most practical navigation and mapping work.
Where haversine is not enough:
- High-precision surveying, cadastral mapping, or scientific geodesy require ellipsoidal methods such as Vincenty’s formulae or Karney’s geodesic algorithm (used in GeographicLib).
- For routing over roads, the dominant error usually comes from the fact that roads do not follow great-circle arcs, not from the spherical assumption.
So the key question is not “Is haversine perfect?” but “Is spherical great-circle distance the right approximation for this system?”
Numerical stability notes
The atan2 form is preferred:
rather than:
\[c = 2\arcsin(\sqrt{a})\]Both are mathematically equivalent, but atan2 tends to behave better near the
edges of floating-point precision.
Also note that due to rounding, a can occasionally drift slightly above 1 for
nearly antipodal points. In defensive production code, it is reasonable to
clamp:
const safeA = Math.min(1, Math.max(0, a));
before evaluating the final angle.

Engraved title page of The Advancement and Proficience of LearningPublic Domain
Engraved title page of Bacon’s Novum OrganumPublic Domain



