Computational formulae

On computer systems with low floating point precision, the spherical law of cosines formula can have large rounding errors if the distance is small (if the two points are a kilometer apart on the surface of the Earth, the cosine of the central angle is near 0.99999999). For modern 64-bit floating-point numbers, the spherical law of cosines formula, given above, does not have serious rounding errors for distances larger than a few meters on the surface of the Earth.^[3] The haversine formula is numerically better-conditioned for small distances by using the chord-length relation:^[4]

${\displaystyle {\begin{aligned}\Delta \sigma &=\operatorname {archav} \left(\operatorname {hav} \left(\Delta \phi \right)+\left(1-\operatorname {hav} (\Delta \phi )-\operatorname {hav} (\phi _{1}+\phi _{2})\right)\operatorname {hav} \left(\Delta \lambda \right)\right).\end{aligned}}}$

Historically, the use of this formula was simplified by the availability of tables for the haversine function: hav(θ) = sin²(θ/2).

The following shows the equivalent formula expressing the chord length explicitly:

${\displaystyle {\begin{aligned}\Delta \sigma &=2\arcsin {\frac {\Delta \sigma _{\text{c}}}{2}},\\\Delta \sigma _{\text{c}}&=2{\sqrt {\sin ^{2}{\frac {\Delta \phi }{2}}+\left(\cos ^{2}\phi _{\text{m}}-\sin ^{2}{\frac {\Delta \phi }{2}}\right)\sin ^{2}{\frac {\Delta \lambda }{2}}}},\end{aligned}}}$

where ${\displaystyle \phi _{\text{m}}={\frac {\phi _{1}+\phi _{2}}{2}}}$.

Although this formula is accurate for most distances on a sphere, it too suffers from rounding errors for the special (and somewhat unusual) case of antipodal points. A formula that is accurate for all distances is the following special case of the Vincenty formula for an ellipsoid with equal major and minor axes:^[5]

${\displaystyle \Delta \sigma =\arctan {\frac {\sqrt {\left(\cos \phi _{2}\sin \Delta \lambda \right)^{2}+\left(\cos \phi _{1}\sin \phi _{2}-\sin \phi _{1}\cos \phi _{2}\cos \Delta \lambda \right)^{2}}}{\sin \phi _{1}\sin \phi _{2}+\cos \phi _{1}\cos \phi _{2}\cos \Delta \lambda }}.}$

Here the quadrant for ${\displaystyle \Delta \sigma }$ should be governed by the signs of the numerator and denominator of the right hand side, e.g., using the atan2 function.

Vector version

Another representation of similar formulas, but using normal vectors instead of latitude and longitude to describe the positions, is found by means of 3D vector algebra, using the dot product, cross product, or a combination:^[6]

${\displaystyle {\begin{aligned}\Delta \sigma &=\arccos \left(\mathbf {n} _{1}\cdot \mathbf {n} _{2}\right)\\&=\arcsin \left|\mathbf {n} _{1}\times \mathbf {n} _{2}\right|\\&=\arctan {\frac {\left|\mathbf {n} _{1}\times \mathbf {n} _{2}\right|}{\mathbf {n} _{1}\cdot \mathbf {n} _{2}}}\\\end{aligned}}}$

where ${\displaystyle \mathbf {n} _{1}}$ and ${\displaystyle \mathbf {n} _{2}}$ are the normals to the sphere at the two positions 1 and 2. Similarly to the equations above based on latitude and longitude, the expression based on arctan is the only one that is well-conditioned for all angles. The expression based on arctan requires the magnitude of the cross product over the dot product.

From chord length

A line through three-dimensional space between points of interest on a spherical Earth is the chord of the great circle between the points. The central angle between the two points can be determined from the chord length. The great circle distance is proportional to the central angle.

The great circle chord length, ${\displaystyle \Delta \sigma _{\text{c}}\,\!}$, may be calculated as follows for the corresponding unit sphere, by means of Cartesian subtraction:

${\displaystyle {\begin{aligned}\Delta {X}&=\cos \phi _{2}\cos \lambda _{2}-\cos \phi _{1}\cos \lambda _{1};\\\Delta {Y}&=\cos \phi _{2}\sin \lambda _{2}-\cos \phi _{1}\sin \lambda _{1};\\\Delta {Z}&=\sin \phi _{2}-\sin \phi _{1};\\\Delta \sigma _{\text{c}}&={\sqrt {(\Delta {X})^{2}+(\Delta {Y})^{2}+(\Delta {Z})^{2}}}.\end{aligned}}}$