The haversine formula is a mathematical equation used to compute the great-circle distance between two points on the surface of a sphere, such as Earth, given their latitudes and longitudes.[1] It relies on the haversine function, defined as \operatorname{hav}(\theta) = \sin^2(\theta/2) = \frac{1 - \cos \theta}{2}, which is half the versine of \theta and offers numerical stability for small angular differences by avoiding catastrophic cancellation in trigonometric computations.[2] This formula is particularly valuable in navigation, geodesy, and geographic information systems (GIS) for accurately modeling distances on spherical approximations of planetary bodies, assuming a uniform radius like Earth's mean value of approximately 6,371 km.[3]
The haversine function originated in the context of spherical trigonometry and logarithmic tables for astronomical calculations. It first appeared in printed tables of logarithmic versines compiled by Spanish mathematician and astronomer José de Mendoza y Ríos in 1801, facilitating computations for lunar distance methods in determining longitude at sea. The term "haversine," a contraction of "half-versed-sine," was coined in 1835 by British mathematician and clergyman James Inman in the third edition of his textbook Navigation and Nautical Astronomy, where he introduced it to simplify nautical calculations involving small angles.[4]
In its standard form for distance calculation, the haversine formula is expressed as a = \operatorname{hav}(\phi_2 - \phi_1) + \cos \phi_1 \cos \phi_2 \operatorname{hav}(\lambda_2 - \lambda_1), where \phi_1, \phi_2 are the latitudes and \lambda_1, \lambda_2 are the longitudes (in radians), followed by the central angle \Delta \sigma = 2 \operatorname{asin}(\sqrt{a}) and distance d = R \Delta \sigma, with R as the sphere's radius.[1] This arrangement derives from the spherical law of cosines but uses haversines to enhance precision for short distances, as highlighted in navigational literature.[3] Though superseded by more precise ellipsoidal models like Vincenty's formula for high-accuracy geodesy, the haversine formula remains a foundational tool due to its simplicity and effectiveness for spherical approximations in applications ranging from aviation routing to mobile app geolocation.[3]
Fundamentals
Haversine Function
The haversine function, denoted \operatorname{hav}(\theta), is defined mathematically as
\operatorname{hav}(\theta) = \sin^2\left(\frac{\theta}{2}\right),
where \theta is an angular measure typically expressed in radians.[2] This formulation arises from half-angle identities in trigonometry and provides a convenient way to express certain spherical trigonometric relations.
The haversine is closely related to the versine function, defined as \operatorname{versin}(\theta) = 1 - \cos(\theta), such that
\operatorname{hav}(\theta) = \frac{\operatorname{versin}(\theta)}{2}.
[2] This equivalence highlights its connection to cosine-based computations, but the sine-squared form is often preferred for its numerical stability.
For \theta \in [0, \pi], the haversine function is continuous, increasing, and bounded between 0 and 1, with \operatorname{hav}(0) = 0 and \operatorname{hav}(\pi) = 1.[2] It proves especially valuable in avoiding catastrophic cancellation errors during small-angle approximations; for instance, when \theta is small, \cos(\theta) \approx 1 - \frac{\theta^2}{2}, but direct subtraction of nearly equal values leads to precision loss, whereas \sin^2(\theta/2) \approx (\theta/2)^2 remains accurate using standard floating-point arithmetic.
The haversine function was first utilized in navigational tables by Spanish astronomer José de Mendoza y Ríos in 1801, as documented by mathematician and historian Florian Cajori. The term "haversine," a contraction of "half-versed-sine," was introduced in 1835 by English mathematician James Inman in his textbook Navigation and Nautical Astronomy to aid spherical calculations in astronomy and navigation. This function underpins the haversine formula for great-circle distances between points on a sphere.
Spherical Geometry Basics
Spherical coordinates provide a system for locating points on the surface of a sphere using angular measurements relative to a reference plane and axis. Latitude, denoted as φ, measures the angle north or south of the equator, ranging from -90° at the South Pole to +90° at the North Pole, while longitude, denoted as λ, measures the angle east or west from the Prime Meridian, ranging from -180° to +180°.[5] These coordinates are defined with respect to the center of the sphere, treating latitude and longitude as angles on a sphere of radius R.[6]
The great-circle distance represents the shortest path between two points on a sphere, corresponding to the arc of a great circle that passes through both points and the sphere's center. This distance is given by d = R \Delta \sigma, where \Delta \sigma is the central angle subtended by the two points at the sphere's center, measured in radians.[7][8] Great circles are the largest possible circles on the sphere, with the same radius as the sphere itself, and they form the basis for navigation and geodesic calculations.[7]
In applications involving Earth, the planet is often modeled as a perfect sphere with a mean radius of approximately 6371 km to simplify computations, though this ignores the ellipsoidal shape and variations due to rotation and topography.[9] This spherical approximation is sufficient for many purposes but introduces minor errors compared to more precise ellipsoidal models.[9]
Points on the sphere can be represented as vectors in three-dimensional Cartesian coordinates for facilitating angular separation calculations. For a point at latitude φ and longitude λ on a unit sphere (where R = 1), the position vector is (\cos \phi \cos \lambda, \cos \phi \sin \lambda, \sin \phi).[10] Scaling by the sphere's radius R yields the full coordinates: (R \cos \phi \cos \lambda, R \cos \phi \sin \lambda, R \sin \phi).[10] These vector forms enable dot product computations to determine the central angle \Delta \sigma via \cos \Delta \sigma = \mathbf{p_1} \cdot \mathbf{p_2} / R^2, providing a foundation for distance-related trigonometric functions.[10]
Distance Calculation
The Haversine formula computes the great-circle distance, the shortest path along the surface of a sphere between two points specified by their latitudes and longitudes. It is particularly suited for numerical stability in calculations involving small angular separations, making it valuable in navigation and geospatial applications. The formula assumes a spherical model of the Earth, ignoring its oblateness for simplicity.[11]
The required inputs are the latitudes \phi_1 and \phi_2 of the two points, and the absolute difference in longitudes \Delta\lambda = |\lambda_2 - \lambda_1|, with all angles expressed in radians. Longitudes and latitudes are typically provided in degrees and must be converted to radians by multiplying by \pi / 180. The formula proceeds in steps to avoid direct computation of small differences that could lead to loss of precision:
\begin{align*}
a &= \sin^2\left(\frac{\Delta\phi}{2}\right) + \cos \phi_1 \cdot \cos \phi_2 \cdot \sin^2\left(\frac{\Delta\lambda}{2}\right), \\
c &= 2 \cdot \atan2\left(\sqrt{a}, \sqrt{1 - a}\right), \\
d &= R \cdot c,
\end{align*}
where R is the radius of the sphere, and the terms \sin^2(\theta/2) correspond to the haversine function \hav(\theta).[11]
The output d represents the great-circle distance in the same units as R. For Earth-distance calculations, R is often taken as the mean radius of 6371 km; alternatively, the equatorial radius of 6378 km may be used for greater accuracy near the equator. These values derive from precise geodetic measurements.[12]
Law of Haversines
The law of haversines is a fundamental relation in spherical trigonometry that expresses any side or angle of a spherical triangle in terms of the haversines of the other elements. It applies to spherical triangles formed by the intersections of great circles on a sphere, where the sides represent angular distances and the angles are dihedral angles between the planes of the great circles.[13]
A common form of the law for the side c opposite angle C in a spherical triangle with sides a, b and angles A, B is given by
\hav c = \hav(a - b) + \sin a \sin b \, \hav C
where \hav \theta = \frac{1 - \cos \theta}{2} = \sin^2 \left( \frac{\theta}{2} \right).[14] This form is particularly useful for navigation computations involving known sides and the included angle.
Analogous expressions exist for the other sides and for the angles, enabling symmetric treatment of all elements in the triangle. These formulas facilitate solving oblique spherical triangles without directly computing cosines, which is advantageous in tabular computations.
Compared to the spherical law of cosines, the law of haversines offers improved numerical stability, particularly for small angular distances where cosine values are close to unity, leading to subtractive cancellation errors; the haversine function, being quadratic in small angles, preserves precision in such cases.[15] This property made it valuable in pre-digital navigation and persists in modern computational geodesy for accurate triangle solutions.
The haversine distance formula between two points on a sphere arises as a special case of the law of haversines when the spherical triangle degenerates into a configuration defined by the points' latitudes and longitude difference. This extends the basic two-point distance calculation to more complex spherical geometries involving multiple great-circle arcs.
Derivation
Geometric Proof
To derive the Haversine formula geometrically, represent the positions of two points on a unit sphere as unit vectors in three-dimensional Cartesian coordinates, assuming a spherical model of Earth. Let the points P_1 and P_2 have latitudes \phi_1, \phi_2 and longitudes \lambda_1, \lambda_2, where latitudes range from -\pi/2 to \pi/2 and longitudes from -\pi to \pi. The corresponding unit vectors are given by:
\begin{align*}
\mathbf{P_1} &= \begin{pmatrix} \cos\phi_1 \cos\lambda_1 \\ \cos\phi_1 \sin\lambda_1 \\ \sin\phi_1 \end{pmatrix}, \\
\mathbf{P_2} &= \begin{pmatrix} \cos\phi_2 \cos\lambda_2 \\ \cos\phi_2 \sin\lambda_2 \\ \sin\phi_2 \end{pmatrix}.
\end{align*}
[1]
The central angle \Delta\sigma subtended by the arc between P_1 and P_2 at the sphere's center equals the angle between the vectors \mathbf{P_1} and \mathbf{P_2}, which is found using the dot product formula for unit vectors:
\cos \Delta\sigma = \mathbf{P_1} \cdot \mathbf{P_2}.
[1]
Expanding the dot product explicitly produces the spherical law of cosines in terms of latitude and longitude difference \Delta\lambda = \lambda_2 - \lambda_1:
\cos \Delta\sigma = \sin \phi_1 \sin \phi_2 + \cos \phi_1 \cos \phi_2 \cos \Delta\lambda.
[1]
To arrive at the Haversine formula, apply the half-angle identity for the haversine function, defined as \hav \theta = \sin^2(\theta/2) = (1 - \cos \theta)/2:
$1 - \cos \Delta\sigma = 2 \sin^2 \left( \frac{\Delta\sigma}{2} \right),
so
\hav(\Delta\sigma) = \frac{1 - \cos \Delta\sigma}{2}.
Substituting the expanded cosine expression and applying trigonometric half-angle identities \sin^2(\alpha/2) = (1 - \cos \alpha)/2 to the latitude difference \Delta\phi = \phi_2 - \phi_1 and longitude difference terms yields:
\hav(\Delta\sigma) = \sin^2 \left( \frac{\Delta\phi}{2} \right) + \cos \phi_1 \cos \phi_2 \sin^2 \left( \frac{\Delta\lambda}{2} \right).
This equation gives the haversine of the central angle directly in terms of the input coordinates, avoiding issues with cosine near unity for small angles.[1][16]
Geometrically, this derivation corresponds to the spherical triangle formed by the two points P_1, P_2, and the north pole N, where the sides from N to P_1 and N to P_2 are the co-latitudes \pi/2 - \phi_1 and \pi/2 - \phi_2, the angle at N is \Delta\lambda, and the side P_1P_2 has length \Delta\sigma. The vector approach embeds this triangle in the plane through the sphere's center, using the dot product to capture the angular separation.[1]
Trigonometric Derivation
The trigonometric derivation of the haversine formula relies on the spherical law of cosines, a fundamental relation in spherical trigonometry for any spherical triangle with sides a, b, c (in angular measure) and opposite angles A, B, C:
\cos c = \cos a \cos b + \sin a \sin b \cos C
This equation, established in classical treatments of spherical geometry, allows computation of one side given the other two sides and the included angle.
To derive the great-circle distance between two points on a sphere specified by latitudes \phi_1 and \phi_2 (in radians) and longitude difference \Delta\lambda, construct a spherical triangle with vertices at the two points and the north pole. The sides adjacent to the north pole are the co-latitudes: a = \frac{\pi}{2} - \phi_1 and b = \frac{\pi}{2} - \phi_2. The angle C at the north pole equals \Delta\lambda.
Substituting these into the spherical law of cosines yields:
\cos c = \sin\phi_1 \sin\phi_2 + \cos\phi_1 \cos\phi_2 \cos\Delta\lambda
Here, c represents the central angle subtended by the great-circle arc between the points; multiplying by the sphere's radius R gives the distance d = R c. This form directly follows from the trigonometric substitutions \cos a = \sin\phi_1, \sin a = \cos\phi_1, and analogously for b.
To express this in haversine form, apply the haversine identity \hav \theta = \frac{1 - \cos \theta}{2} = \sin^2 \frac{\theta}{2}, which rearranges the law of cosines into a form suited for navigational tables and avoiding subtraction of near-equal quantities for small angles. The general identity transforming the spherical law of cosines to haversines is the law of haversines:
\hav c = \hav(a - b) + \sin a \sin b \hav C
[13]
Substituting a = \frac{\pi}{2} - \phi_1, b = \frac{\pi}{2} - \phi_2, and C = \Delta\lambda simplifies to the haversine formula, since \sin a = \cos \phi_1, \sin b = \cos \phi_2, \hav(a - b) = \hav(\phi_2 - \phi_1), and \hav C = \hav(\Delta\lambda):
\hav c = \hav(\phi_2 - \phi_1) + \cos\phi_1 \cos\phi_2 \hav(\Delta\lambda)
This equation, a special case of the broader law of haversines in spherical trigonometry, facilitates accurate distance computation, particularly for small separations where cosine-based forms suffer numerical instability.
Applications and Examples
Navigation and Geodesy
The Haversine formula has been historically employed in celestial navigation to compute great-circle distances between observed celestial bodies and a vessel's position, facilitating sight reductions that integrate with dead reckoning techniques for position fixes prior to the advent of GPS.[17][18] In dead reckoning, navigators used the formula to estimate distances traveled based on course, speed, and time, often in conjunction with haversine tables for logarithmic computations to minimize errors over long voyages.[17] Its adoption extended to aviation route planning in the early 20th century, where it enabled pilots to calculate efficient great-circle paths for transoceanic flights, reducing fuel consumption and flight time on spherical Earth models.[19]
In modern applications, the Haversine formula serves as a foundational method for distance calculations between positions obtained from GPS systems, such as in smartphones, by determining great-circle distances between latitude-longitude coordinates.[11] Many mapping applications use the Haversine formula for straight-line distance calculations and proximity queries, often in conjunction with APIs like Google Maps for displaying results on spherical Earth approximations.[11]
Within geodesy, the formula approximates distances on a spherical model of Earth, treating the planet as a perfect sphere with a mean radius of approximately 6371 km, which simplifies computations for global-scale analyses but introduces errors when compared to more precise ellipsoidal representations.[11] For high-precision requirements, such as those using the WGS84 ellipsoid standard in satellite geodesy, the Haversine formula's limitations necessitate alternatives like Vincenty's formula, which accounts for Earth's oblate shape and yields geodesic distances accurate to within millimeters over long baselines.[20]
Error analysis reveals that the Haversine formula achieves an accuracy of about 0.3% to 0.5% relative to true geodesic distances on Earth, making it suitable for most navigational purposes but less ideal for sub-millimeter geodetic surveys.[21][20] Its numerical stability is particularly advantageous for small distances, where it avoids subtraction cancellation errors inherent in other spherical trigonometric approaches, ensuring reliable results even for separations under 1 km.[11]
Computational Implementation
The computational implementation of the Haversine formula involves converting latitude and longitude coordinates to radians, computing the haversine of the central angle, and applying the inverse haversine to obtain the great-circle distance.[11]
A standard pseudocode outline for calculating the distance d between two points with latitudes \phi_1, \phi_2 and longitudes \lambda_1, \lambda_2 (in degrees) on a sphere of radius R is as follows:
function haversine_distance(φ1_deg, λ1_deg, φ2_deg, λ2_deg, [R](/page/R)):
# Convert degrees to radians
φ1 = φ1_deg * π / 180
λ1 = λ1_deg * π / 180
φ2 = φ2_deg * π / 180
λ2 = λ2_deg * π / 180
# Differences in radians
Δφ = φ2 - φ1
Δλ = λ2 - λ1
# Haversine computation
a = sin²(Δφ / 2) + cos(φ1) * cos(φ2) * sin²(Δλ / 2)
# Inverse haversine using [atan2](/page/Atan2) for [angular distance](/page/Angular_distance) c
c = 2 * [atan2](/page/Atan2)(√a, √(1 - a))
# [Distance](/page/Distance)
d = R * c
return d
function haversine_distance(φ1_deg, λ1_deg, φ2_deg, λ2_deg, [R](/page/R)):
# Convert degrees to radians
φ1 = φ1_deg * π / 180
λ1 = λ1_deg * π / 180
φ2 = φ2_deg * π / 180
λ2 = λ2_deg * π / 180
# Differences in radians
Δφ = φ2 - φ1
Δλ = λ2 - λ1
# Haversine computation
a = sin²(Δφ / 2) + cos(φ1) * cos(φ2) * sin²(Δλ / 2)
# Inverse haversine using [atan2](/page/Atan2) for [angular distance](/page/Angular_distance) c
c = 2 * [atan2](/page/Atan2)(√a, √(1 - a))
# [Distance](/page/Distance)
d = R * c
return d
This implementation relies on the two-argument arctangent function (atan2) to ensure the correct quadrant and avoid division by zero in the inverse haversine step.[11]
The Haversine formula offers improved numerical stability compared to the spherical law of cosines, particularly for points in close proximity where the central angle is small. In the law of cosines, computing \cos c = \cos a \cos b + \sin a \sin b \cos C involves subtracting two nearly equal large values (close to 1) when distances are short, leading to catastrophic cancellation and loss of precision in floating-point arithmetic. By contrast, the Haversine approach uses half-angle sine terms, which remain well-conditioned even as \Delta\phi and \Delta\lambda approach zero, preserving accuracy down to sub-millimeter scales on Earth's surface.[11][22]
In programming practice, implementations should employ double-precision floating-point arithmetic (64-bit) to minimize rounding errors, as single-precision (32-bit) may introduce inaccuracies exceeding 1 meter for global distances. Special handling is required for antipodal points, where the central angle \Delta\sigma = \pi; here, \sin(\Delta\lambda / 2) = 0 and \sin(\Delta\phi / 2) = 0, causing the argument to atan2 to be (0, 0) and yielding NaN—in such cases, the distance can be set directly to \pi R. Additionally, input validation for latitude ([-90^\circ, 90^\circ]) and longitude ([-180^\circ, 180^\circ]) ranges is recommended to prevent extraneous results.[11]
The Haversine formula is supported in various programming libraries and geospatial tools. In Python, the haversine package provides a straightforward function for distance computation using the formula.[23] For JavaScript, the geolib library includes Haversine-based distance calculations suitable for web applications.
Numerical Example
To illustrate the application of the Haversine formula, consider calculating the great-circle distance between New York City at 40.7128° N, 74.0060° W and London at 51.5074° N, 0.1278° W.[24][25]
The Haversine formula for the central angle c between two points with latitudes \phi_1, \phi_2 and longitudes \lambda_1, \lambda_2 (in radians) is given by
a = \sin^2\left(\frac{\Delta\phi}{2}\right) + \cos\phi_1 \cos\phi_2 \sin^2\left(\frac{\Delta\lambda}{2}\right),
c = 2 \arcsin(\sqrt{a}),
where \Delta\phi = \phi_2 - \phi_1 and \Delta\lambda = \lambda_2 - \lambda_1, followed by the distance d = R c using Earth's mean radius R = 6371 km.[11]
Convert the coordinates to radians: \phi_1 \approx 0.7111, \phi_2 \approx 0.8988, \lambda_1 \approx -1.2922, \lambda_2 \approx -0.0022. Then, \Delta\phi \approx 0.1877 and \Delta\lambda \approx 1.2900.
Compute \sin^2(\Delta\phi/2) \approx 0.0088, \cos\phi_1 \approx 0.7580, \cos\phi_2 \approx 0.6216, and \sin^2(\Delta\lambda/2) \approx 0.3625. Thus, a \approx 0.0088 + (0.7580 \times 0.6216 \times 0.3625) \approx 0.0088 + 0.1709 = 0.1797.
Next, \sqrt{a} \approx 0.4238, so c = 2 \arcsin(0.4238) \approx 2 \times 0.437 = 0.874 radians. The distance is d \approx 6371 \times 0.874 \approx 5570 km.[11]
This result aligns closely with the established great-circle distance of approximately 5569 km between these points, with the minor discrepancy attributable to the formula's assumption of a perfect sphere rather than Earth's oblate spheroid shape.[11] For greater precision in specific regions, R can be adjusted: use the equatorial radius of 6378 km for paths near the equator or the polar radius of 6357 km for polar routes.[26]