Distances in the Sky

1. Relative Inaccuracy ... 2. Polar Coordinates ... 2.1. Formulas for the Distance Between Polar Coordinates ... 3. Cartesian Coordinates

\( \DeclareMathOperator{\del}{∆\!} \)

This page explains how to calculate distances in the sky. There are many different formulas for that. Those formulas all yield mathematically equivalent results, but for use in numerical calculations (with calculating machine or computer program) some are better than others. So, we first take a look at what makes such a formula good or bad for numerical use.

As an example, we'll calculate the distance between a point with latitude 50.850° and longitude 4.350° and a point with latitude 52.383° and longitude 4.900°. (On Earth, these geographical coordinates belong to Brussels and Amsterdam.)

In addition, we calculate the distance between point A and point C with latitude 50.850001° and longitude 4.900001°, very close to point A.

If you want to use a formula to calculate a distance numerically, then you need to take into account the limited precision of your calculating machine or computer program. There are many formulas that, mathematically speaking, all calculate exactly the same distance. Some of the above equations then yield results for some distances that are much less accurate than you might expect, while the other equations perform fine. Why is that?

Suppose that you have a formula that yields a function of the angular distance \( d \):

\begin{equation} f(d) = q \end{equation}

where \( f \) is a function of only \( d \), and \( q \) is a (possibly complicated) function of the coordinates of the two points. If you know the inverse function \( f^{\text{inv}} \) of \( f \), for which

\begin{equation} f^{\text{inv}}(f(x)) = x \end{equation}

then

\begin{equation} d = f^{\text{inv}}(q) \end{equation}

Now suppose that you shift one of the two points by a small amount. Then \( q \) changes to \( q + \del q \), and \( d \) changes to \( d + \del d \), with

\begin{equation} f'(d) \del d ≈ \del q \end{equation}

where \( f'(d) \) is the derivative of \( f(d) \) with respect to \( d \). The relative difference is

\begin{equation} \frac{\del q}{q} = \frac{f'(d)\del d}{f(d)} = \frac{f'(d)d}{f(d)} \frac{\del d}{d} = N \frac{\del d}{d} \end{equation}

What if \( N = f'(d)d/f(d) \) is much smaller than 1, or perhaps even equal to 0? Then, relatively speaking, \( q \) changes hardly at all when \( d \) does change considerably.

A calculating machine (or computer program) usually calculates with a limited relative accuracy. If two numbers have a relative difference that is smaller than a certain limiting value \( p \) (which depends on the calculating machine), then that calculating machine cannot make that distinction and regards both numbers as exactly the same. If \( \del q/q = p \), then \( \del q \) is the smallest absolute difference that the calculating machine can handle for numbers near \( q \). The corresponding relative difference in \( d \) is then \( p/N \), which is much greater than \( p \) if \( N \) is much smaller than 1, so then \( d \) is known with much less accuracy than you would expect based on \( p \).

To get the greatest accuracy for \( d \), we need to use an equation for which \( |N| \) is not much smaller than 1 for all distances (between 0° and 180°, an angular distance cannot be outside that range). How do the first functions that come to mind, \( \sin, \cos, \tan \) fare?

\({f(d)}\) | \({N}\) | \({\{d | N = 0\}}\) |
---|---|---|

\({\sin(d)}\) | \({d\cos(d)/\sin(d)}\) | 90° |

\({\cos(d)}\) | \({d\sin(d)/\cos(d)}\) | 0°, 180° |

\({\tan(d)}\) | \({d/(\sin(d)\cos(d))}\) | |

The third column says for which distances \( d \) from 0° through 180° the value of \( N \) is equal to 0.

So, a sine yields relatively inaccuracte results for distances near 90°; a cosine for distances near 0° and 180°, and a tangent is relatively accurate for all distances.

\( f^n(d) \) for \( n \gt 0 \) yields the same roots for \( N \) (values of \( d \) for which \( N \) is zero) as \( f(d) \) does. For example, \( \sin^2(d) \) yields relatively inaccurate results for distances near 90°, just like \( \sin(d) \) does.

\( f(kd) \) for constant \( k \gt 0 \) yields roots for the same distances as \( f(d) \) does, but divided by \( k \). For example, \( \sin\left( \frac{1}{2}d \right) \) yields relatively inaccurate results for distances near 180°, twice as large as the 90° of \( \sin(d) \).

So, we can move the root of \( N \) of a sine outside of the interval from 0° to 180° if we can make \( k \) less than 1/2, but that doesn't work for a cosine, because the root at 0° remains there. It is better to find a formula for \( \tan(d) \), because \( N \) for that function has no roots, so \( d \) can be accurately found from \( q \) for all distances.

A disadvantage of \( \tan(d) = q \) itself is that its inverse function \( d = \arctan(q) \) for negative \( q \) yields values \( −90° \lt d \lt 0° \) on most calculating machines. You need to add 180° by hand (which is allowed, because \( \tan(d) = \tan(d + 180°) \), or you need to use the "atan2" function with two arguments (the numerator and denominator of \( q \)), instead of the "tan" function.

In that regard, \( \tan\left( \frac{1}{2}d \right) = q \) is better, because then \( q \ge 0 \) for all \( 0° \le d \le 180° \), so there is no longer a problem for negative \( q \).

Polar coordinates are coordinates of which one is similar to the geographical longitude and the other to the geographical latitude. (If the coordinates refer to a location in space, then the linear distance from the origin of the coordinate system should also be specified, but that distance is not relevant here because we discuss positions in the sky.) In polar coordinates, the latitude measures the distance from the equator, and the longitude measures the distance along the equator. The poles are then at latitudes +90° and −90°, and the equator is at latitude 0°.

Examples of polar coordinates are: longitude and latitude (geographical or equatorial or ecliptic or galactic or celestial body + "graphic"), or right ascension and declination, or azimuth and height.

If you have polar coordinates \( l_1 \) and \( l_2 \) (the longitudes) and \( b_1 \) and \( b_2 \) (the latitudes) of two directions in the sky, then there are many different formulas with which you can calculate the angular distance between those two points on the sky globe. Here is a collection of such formulas:

\begin{align} \cos(d) & = \sin(b_1)\sin(b_2) + \cos(b_1)\cos(b_2)\cos(l_2 − l_1) \label{eq:cos} \\ \cos^2\left( \frac{1}{2}d \right) & = \sin^2\left( \frac{1}{2} (b_2 + b_1) \right) + \cos(b_1)\cos(b_2)\cos^2\left( \frac{1}{2} (l_2 − l_1) \right) \label{eq:cossqh} \\ \sin^2(d) & = \left(\cos(b_2) \sin(l_2 − l_1) \right)^2 + \left( \cos(b_1)\sin(b_2) − \sin(b_1)\cos(b_2)\cos(l_2 − l_1) \right)^2 \label{eq:sinsq} \\ \sin^2\left( \frac{1}{2}d \right) & = \sin^2\left( \frac{1}{2} (b_2 − b_1) \right) + \cos(b_1)\cos(b_2)\sin^2\left( \frac{1}{2} (l_2 − l_1) \right) \label{eq:sinsqh} \\ \tan(d) & = \frac{\sqrt{\left( \cos(b_2) \sin(l_2 − l_1) \right)^2 + \left( \cos(b_1)\sin(b_2) − \sin(b_1)\cos(b_2)\cos(l_2 − l_1) \right)^2}}{\sin(b_1)\sin(b_2) + \cos(b_1)\cos(b_2)\cos(l_2 − l_1)} \label{eq:tan} \\ \tan^2\left( \frac{1}{2}d \right) & = \frac{\sin^2\left( \frac{1}{2} (b_2 − b_1) \right) + \cos(b_1)\cos(b_2)\sin^2\left( \frac{1}{2} (l_2 − l_1) \right)}{\sin^2\left( \frac{1}{2} (b_2 + b_1) \right) + \cos(b_1)\cos(b_2)\cos^2\left( \frac{1}{2} (l_2 − l_1) \right)} \label{eq:tansqh} \end{align}

The transformation into \( d \) then proceeds as follows:

\({ q }\) | \({ d }\) |
---|---|

\({ \cos(d) }\) | \({ \arccos(q) }\) |

\({ \cos^2(\frac{1}{2}d) }\) | \({ 2\arccos(\sqrt{q}) }\) |

\({ \sin^2(d) }\) | \({ \arcsin(\sqrt{q}) }\) |

\({ \sin^2(\frac{1}{2}d) }\) | \({ 2\arcsin(\sqrt{q}) }\) |

\({ \tan(d) }\) | \({ \arctan(d) }\) |

\({ \tan^2(\frac{1}{2}d) }\) | \({ 2\arctan(\sqrt{q}) }\) |

Equation \eqref{eq:cos} is nowadays the best known one, and follows from the spherical law of cosines. Equation \eqref{eq:sinsqh} is the haversine formula and was much used before electronic calculating machines existed. Equations \eqref{eq:sinsq} and \eqref{eq:tan} are special cases of the Vincenty formula from geodesy. Note that the denominator of equation \eqref{eq:tan} is equal to the right-hand side of the cosine formula (Equation \eqref{eq:cos}). I derived Equations \eqref{eq:cossqh} and \eqref{eq:tansqh}. Note that the denominator of Equation \eqref{eq:tansqh} is equal to the right-hand side of the haversine formula (Equation \eqref{eq:sinsqh}).

A disadvantage of Equation \eqref{eq:tan} is that its right-hand side is negative for \( d \gt 90° \) and that the \( \arctan \) function in most calculating machines and computer programs produces a result between −90° and 0° for a negative argument, instead of the result between 90° and 180° that we need. You then have to remember that \( \tan(d) = \tan(d + 180°) \), so you can add 180° to the result. Or you can use the "atan2" function (arctan with two arguments), if it is available, with the numerator and denominator of Equation \eqref{eq:tan} for arguments. The order of the arguments depends on the calculating machine. For \( q = \frac{0}{1} \) the result should be \( d = 0° \), and for \( q = \frac{1}{0} \) the result should be \( d = 90° \).

Mathematically, the above equations are all equally good, but for numerical use they aren't. I've sorted the above list of equations in order of increasing suitability for use in numerical calculations (with calculating machine or computer program), so I think that Equation \eqref{eq:tansqh} is the most suitable for numerical calculations, and Equation \eqref{eq:cos} is the least suitable. I based my judgement mostly on the accuracy of the results, and also on the convenience of calculation.

Equation \eqref{eq:tansqh} mentions both \( \sin^2\left( \frac{1}{2} (l_2 − l_1) \right) \) and \( \cos^2\left( \frac{1}{2} (l_2 − l_1) \right) \), which sum to exactly 1, so you might be tempted to only calculate one of them explicitly and then find the other one by subtracting the first one from 1, but if you do that then the results will in some cases be a lot less accurate.

For points A and B, and using a computer system that works with about 16 significant figures, Equation \eqref{eq:cos} yields:

\begin{align*} q & = \sin(50.850°) \sin(52.383°) + \cos(50.850°) \cos(52.383°) \cos(4.350° − 4.900°) \\ & = 0.775495743172234×0.792108574769034 + 0.631352795449378×0.610380214110329×0.999953926969049 \\ & = 0.999624327383860 \\ d & = \arccos(0.999624327383860) = 1.57056529603582° \end{align*}

Equation \eqref{eq:cossqh} yields:

\begin{align*} q & = \sin^2\left( \frac{52.383° + 50.850°}{2} \right) + \cos(50.850°)\cos(52.383°)\cos^2\left( \frac{4.350° − 4.900°}{2} \right) \\ & = 0.614455786699033 + 0.631352795449378×0.610380214110329×0.999976963484524 \\ & = 0.999812163691930 \\ d & = 2\arccos{\sqrt{0.999812163691930}} = 1.57056529603563° \end{align*}

Equation \eqref{eq:sinsq} yields:

\begin{align*} q & = \sin^2\left( \frac{52.383° − 50.850°}{2} \right) + \cos(50.850°) \cos(52.383°) \sin^2\left( \frac{4.900° − 4.350°}{2} \right) \\ & = (0.610380214110329×−0.00959916346240015)^2 + (0.631352795449378×0.792108574769034 − 0.775495743172234×0.610380214110329×0.999953926969049)^2 \\ & = 0.000751204102365001 \\ d & = \arcsin(\sqrt{0.000751204102365001}) = 1.57056529603550° \end{align*}

Equation \eqref{eq:sinsqh} yields:

\begin{align*} q & = \sin^2\left( \frac{1}{2} (52.383° − 50.850°) \right) + \cos(50.850°) \cos(52.383°) \sin^2\left( \frac{1}{2} (4.900° − 4.350°) \right) \\ & = 0.000178958835421660 + 0.631352795449378×0.610380214110329×0.0000230365154755148 \\ & = 0.000187836308069881 \\ d & = 2×\arcsin(\sqrt{0.000187836308069881}) = 1.57056529603551° \end{align*}

Equation \eqref{eq:tan} yields:

\begin{align*} q & = \frac{\sqrt{0.000751204102365001}}{0.999624327383860} = 0.0274184032029088 \\ d & = \arctan(0.0274184032029088) = 1.57056529603550° \end{align*}

Equation \eqref{eq:tansqh} yields:

\begin{align*} q & = \frac{0.000187836308069881}{0.999812163691930} = 0.000187871597177086 \\ d & = 2\arctan(\sqrt{0.000187871597177086}) = 1.57056529603551° \end{align*}

Using a program that works with even more significant figures, I found \( d = 1.57056529603550534° \) (with both Equations \eqref{eq:cos} and \eqref{eq:tansqh}). Now we can see how far the results of the various equations were from the correct one, in units of the 14th figure after the decimal marker. For Equation \eqref{eq:cos} that was 31; for Equation \eqref{eq:cossqh} it was 12, and for the other equations it was 1 or less. It seems that Equations \eqref{eq:cos} and \eqref{eq:cossqh} are somewhat less accurate in this case.

For the distance between points A and C we find, for the same formulas,

Eq. | \({ d×10^6 }\) | \({ \del d×10^6 }\) |
---|---|---|

\eqref{eq:cos} | 1.207418270 | 0.025 |

\eqref{eq:cossqh} | 0 | 1.18 |

\eqref{eq:sinsq} | 1.182626877 | 0.0000000018 |

\eqref{eq:sinsqh} | 1.182626880 | 0.0000000054 |

\eqref{eq:tan} | 1.182626877 | 0.0000000018 |

\eqref{eq:tansqh} | 1.182626880 | 0.0000000054 |

\( \del d \) is de deviation with respect to the correct answer \( 1.1826268750×10^{−6}° \). It is clear that Equations \eqref{eq:cos} and \eqref{eq:cossqh} perform much worse here than the other equations do: Equation \eqref{eq:cossqh} cannot distinguish the answer from 0, and the deviation is about 10 million times greater for Equation \eqref{eq:cos} than for the other equations. The equations that perform well yield an answer that is correct to about 14 figures after the decimal marker.

Let's investigate in more detail the relative inaccuracy of Equation \eqref{eq:cos} for small distances. For small \( d \), and if we measure all angles in radians, then, approximately

\begin{align} q & = \cos(d) ≈ 1 − \frac{1}{2}d^2 \\ d & ≈ \sqrt{2(1 − q)} \end{align}

If the smallest relative difference that the calculating machine can handle is \( p \), then the smallest value for \( q \) that can be distinguished from \( q = 1 \) (to which corresponds \( d = 0 \)) is equal to \( q_1 = 1 − p \), for which \( d_1 ≈ \sqrt{2p} \), which is much larger than \( p \).

If your calculating machine or computer program uses uses the "binary64" format from the IEEE 754 standard (also called "floating-point numbers of double precision"), then it calculates with 15 to 17 significant figures. I have such a system, and it yields:

\begin{align*} q_1 & ≈ 0.99999999999999978 \\ d_1 & ≈ 0.000000021 \text{ rad} ≈ 0.0000012° \end{align*}

So the smallest angular distance (other than 0) that you can get out of that formula is then 0.0000012°, which on the surface of the Earth corresponds to 13 cm. Based on \( p ≈ 2×10^{−16} \) you'd expect to be able to handle distances as small as \( 360°×p ≈ 0.0000000000008° \), which on Earth corresponds to 0.000009 mm, much smaller that the 13 cm from before.

Note that the distance between points A and C that followed from Equation \eqref{eq:cos} was equal to \( d_1 \). That distance was therefore at the very limit of what that formula can handle (on that system), and was too far for Equation \eqref{eq:cossqh}.

A calculating machine that uses the "binary32" format (floating-point numbers of single precision) yields

\begin{align*} q_1 & ≈ 0.99999988 \\ d_1 & ≈ 0.00049 \text{ rad} ≈ 0.028° \end{align*}

0.028° on Earth corresponds to about 3 km. Based on \( p ≈ 10^{−7} \) you'd expect to be able to handle distances about as small as \( 360°×p ≈ 0.00004° \), which on Earth corresponds to about 5 m.

If you know the cartesian geocentric coordinates \( x_1, y_1, z_1 \) and \( x_2, y_2, z_2 \) of the two points, then you can calculate their mutual angular distance \( d \) with the formulas mentioned below. In all case we assume that \( x_1^2 + y_1^2 + z_1^2 = 1 \) and \( x_2^2 + y_2^2 + z_2^2 = 1 \). If that isn't true, then divide \( x_1, y_1, z_1 \) by \( \sqrt{x_1^2 + y_1^2 + z_1^2} \), and likewise for \( x_2, y_2, z_2 \).

\begin{align} \cos(d) & = x_1 x_2 + y_1 y_2 + z_1 z_2 \label{eq:rcos} \\ 4\cos^2\left( \frac{1}{2} d \right) & = (x_1 + x_2)^2 + (y_1 + y_2)^2 + (z_1 + z_2)^2 \label{eq:rcossqh} \\ \sin^2(d) & = (y_1 z_2 − z_1 y_2)^2 + (z_1 x_2 − x_1 z_2)^2 + (x_1 y_2 − y_1 x_2)^2 \label{eq:rsinsq} \\ 4\sin^2\left( \frac{1}{2} d \right) & = (x_1 − x_2)^2 + (y_1 − y_2)^2 + (z_1 − z_2)^2 \label{eq:rsinsqh} \\ \tan(d) & = \frac{\sqrt{(y_1 z_2 − z_1 y_2)^2 + (z_1 x_2 − x_1 z_2)^2 + (x_1 y_2 − y_1 x_2)^2}}{x_1 x_2 + y_1 y_2 + z_1 z_2} \label{eq:rtan} \\ \tan^2\left( \frac{1}{2} d \right) & = \frac{(x_1 − x_2)^2 + (y_1 − y_2)^2 + (z_1 − z_2)^2}{(x_1 + x_2)^2 + (y_1 + y_2)^2 + (z_1 + z_2)^2} \label{eq:rtansqh} \end{align}

The transformation into \( d \) then proceeds as follows:

\({ q }\) | \({ d }\) |
---|---|

\({ \cos(d) }\) | \({ \arccos(q) }\) |

\({ 4\cos^2(\frac{1}{2}d) }\) | \({ 2\arccos(\frac{1}{2}\sqrt{q}) }\) |

\({ \sin(d) }\) | \({ \arcsin(q) }\) |

\({ 4\sin^2(\frac{1}{2}d) }\) | \({ 2\arcsin(\frac{1}{2}\sqrt{q}) }\) |

\({ \tan(d) }\) | \({ \arctan(d) }\) |

\({ \tan^2(\frac{1}{2}d) }\) | \({ 2\arctan(\sqrt{q}) }\) |

Equation \eqref{eq:rcos} follows from the inner product of the two vectors. Equation \eqref{eq:rsinsq} follows from the outer product of the two vectors. Equation \eqref{eq:rsinsqh} follows from the length of the straight line connecting the two points (which is the difference between the two vectors). Equation \eqref{eq:rcossqh} follows from the length of the sum of the two vectors. Equation \eqref{eq:rtan} follows from the combination of Equations \eqref{eq:rsinsq} and \eqref{eq:rcos}. Equation \eqref{eq:rtansqh} follows from the combination of Equations \eqref{eq:rsinsqh} and \eqref{eq:rcossqh}. I derived equations \eqref{eq:rcossqh} and \eqref{eq:rtansqh} myself.

Based on the analysis of relative inaccuracy, the formulas that go via \( \tan \) are the best ones. Equation \eqref{eq:rtansqh} is a bit easier to use than Equation \eqref{eq:rtan}, because for Equation \eqref{eq:rtan} you (usually) have to add 180° to the result if \( q \lt 0 \) to get \( 90° \lt d \le 180° \), but for Equation \eqref{eq:rtansqh} that is not necessary.

Using the formulas from section 1.6 of the page about positions in the sky (and \( \del = 1 \)) we find (with a calculating machine that uses about 16 significant figures)

\begin{align*} x_1 & = 0.629045387982967 \\ y_1 & = 0.0539282132014616 \\ z_1 & = 0.775495743172234 \\ x_2 & = 0.608621905592157 \\ y_2 & = 0.0462966717026435 \\ z_2 & = 0.792108574769034 \end{align*}

And then we find from Equation \eqref{eq:rcos}:

\begin{align*} \\ q & = 0.629045387982967 × 0.608621905592157 + 0.0539282132014616 × 0.0462966717026435 + 0.775495743172234 × 0.792108574769034 = 0.999624327383860 \\ d & = \arccos(0.999624327383860) = 1.57056529603582° \end{align*}

Equation \eqref{eq:rcossqh} yields:

\begin{align*} \\ q & = (0.629045387982967 + 0.608621905592157)^2 + (0.0539282132014616 + 0.0462966717026435)^2 + (0.775495743172234 + 0.792108574769034)^2 = 3.99924865476772 \\ d & = 2\arccos\left( \sqrt{\frac{3.99924865476772}{4}} \right) = 1.57056529603656° \end{align*}

Equation \eqref{eq:rsinsq} yields:

\begin{align*} q & = (0.0539282132014616 × 0.792108574769034 − 0.775495743172234 × 0.0462966717026435)^2 \\ & + (0.775495743172234 × 0.608621905592157 − 0.629045387982967 × 0.792108574769034)^2 \\ & + (0.629045387982967 × 0.0462966717026435 − 0.0539282132014616 × 0.608621905592157)^2 \\ & = 0.0000464323440855707 + 0.000691087795493142 + 0.0000136839627862883 \\ & = 0.000751204102365001 \\ d & = \arcsin\left( \sqrt{0.000751204102365001} \right) = 1.57056529603550° \end{align*}

Equation \eqref{eq:rsinsqh} yields:

\begin{align*} \\ q & = (0.629045387982967 − 0.608621905592157)^2 + (0.0539282132014616 − 0.0462966717026435)^2 + (0.775495743172234 − 0.792108574769034)^2 = 0.000751345232279519 \\ d & = 2\arcsin\left( \sqrt{\frac{0.000751345232279519}{4}} \right) = 1.57056529603550° \end{align*}

Equation \eqref{eq:rtan} yields:

\begin{align*} q & = \frac{\sqrt{0.000751204102365001}}{0.999624327383860} = 0.0274184032029088 \\ d & = \arctan(0.0274184032029088) = 1.57056529603550° \end{align*}

Equation \eqref{eq:rtansqh} yields:

\begin{align*} q & = \frac{0.000751345232279519}{3.99924865476772} = 0.000187871597177085 \\ d & = 2\arctan(0.000187871597177085) = 1.57056529603550° \end{align*}

For the distance between points A and C we find, for the same formulas,

Eq. | \({ d×10^6 }\) | \({ |\del d|×10^6 }\) |
---|---|---|

\eqref{eq:rcos} | 1.478779333 | 0.2962 |

\eqref{eq:rcossqh} | 1.707547293 | 0.5249 |

\eqref{eq:rsinsq} | 1.182626877 | 0.0000004 |

\eqref{eq:rsinsqh} | 1.182626877 | 0.0000004 |

\eqref{eq:rtan} | 1.182626877 | 0.0000004 |

\eqref{eq:rtansqh} | 1.182626877 | 0.0000004 |

\( |\del d| \) is the deviation with respect to the correct answer \( 1.182626882738767×10^{−6}° \). Here, like for the polar coordinates, it is clear that the formulas that go via \( \cos \) are a lot less accurate for small distances than the other formulas are.

*//aa.quae.nl/en/reken/afstanden.html;
Last updated: 2017-04-24*