Distance on a sphere — Haversine and friends

Measuring distance between two geographic points is not as simple as applying Pythagoras. Lat/lon coordinates sit on a curved surface, and treating them as flat (x, y) pairs produces errors that grow with distance. At city scale the error is negligible; at country scale it becomes significant. This chapter covers the spectrum from fast-and-approximate to accurate.

Why Pythagoras lies (mostly)

For tiny distances (< ~1 km in mid-latitudes), you can approximate:

$$\Delta x = (\lambda_2 - \lambda_1) \cdot \cos\!\left(\tfrac{\varphi_1 + \varphi_2}{2}\right) \cdot R$$ $$\Delta y = (\varphi_2 - \varphi_1) \cdot R$$ $$d \approx \sqrt{\Delta x^2 + \Delta y^2}$$

This "equirectangular approximation" is fast and accurate to ~0.5% over short distances. Useful in tight inner loops.

The cos(φ) factor in the Δx term corrects for longitude lines converging — without it, a 1° east-west step would appear the same size regardless of latitude, which is only true at the equator.

The Haversine formula

The gold standard for sphere distances:

$$a = \sin^2\!\left(\tfrac{\Delta\varphi}{2}\right) + \cos\varphi_1 \cdot \cos\varphi_2 \cdot \sin^2\!\left(\tfrac{\Delta\lambda}{2}\right)$$ $$c = 2 \cdot \arctan2\!\left(\sqrt{a}, \sqrt{1 - a}\right)$$ $$d = R \cdot c$$

The intermediate variable a is the square of half the chord length between the two points on a unit sphere. The term c is the central angle between them (in radians), and multiplying by R gives the arc length along the surface. The arctan2 form is numerically stable across all distances, including antipodal points — use it instead of 2 * arcsin(sqrt(a)).

In code:

function haversine(lat1, lon1, lat2, lon2) {
  const R = 6371008.8; // meters
  const toRad = d => d * Math.PI / 180;

  const φ1 = toRad(lat1);
  const φ2 = toRad(lat2);
  const Δφ = toRad(lat2 - lat1);
  const Δλ = toRad(lon2 - lon1);

  const a = Math.sin(Δφ/2) ** 2 +
            Math.cos(φ1) * Math.cos(φ2) * Math.sin(Δλ/2) ** 2;
  const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

  return R * c;
}