### How can we find the formula for the great circle distance?

In Lost but lovely: the haversine we explored great circle distances on the Earth and looked at the haversine formula for calculating them: \[\mathop{\mathrm{haversin}} \left(\frac{d}{R}\right) = \mathop{\mathrm{haversin}} (\phi_2-\phi_1) + \cos{\phi_1}\cos{\phi_2} \mathop{\mathrm{haversin}} (\lambda_2-\lambda_1).\]

Recall that in this formula, we are finding the great circle distance \(d\) between two points \(P\) and \(Q\) which have latitudes \(\phi_1\) and \(\phi_2\) and longitudes \(\lambda_1\) and \(\lambda_2\) respectively. The radius of the Earth (or sphere) is \(R\), and the haversine is defined by \[\mathop{\mathrm{haversin}} x=\sin^2\left(\frac{x}{2}\right).\] (On the left hand side of the first equation, \(\frac{d}{R}\) is thought of as an angle in radians.)

So we can rewrite the haversine formula in more familiar notation as \[\sin^2{\left(\frac{d}{2R}\right)} = \sin^2{\left(\frac{\phi_2-\phi_1}{2}\right)} + \cos{\phi_1}\cos{\phi_2} \sin^2{\left(\frac{\lambda_2-\lambda_1}{2}\right)}.\]

In this piece, we explore where this formula comes from, and why navigators might have used the unfamiliar haversine rather than ordinary sines and cosines in their calculations.

#### Preliminaries

We can use any units of distance that we wish when we do our calculations: changing units does not affect the value of \(\frac{d}{R}\) (as long as we use the same units for both \(d\) and \(R\)). Most conveniently, we can choose units so that the radius of the earth is \(R=1\), meaning that we are working with a*unit sphere*. Then the haversine formula becomes \[\begin{equation} \label{eq:haversine} \mathop{\mathrm{haversin}} d = \mathop{\mathrm{haversin}} (\phi_2-\phi_1) + \cos{\phi_1}\cos{\phi_2} \mathop{\mathrm{haversin}} (\lambda_2-\lambda_1) \end{equation}\] or \[\begin{equation} \label{eq:haversinesine} \sin^2{\left(\frac{d}{2}\right)} = \sin^2{\left(\frac{\phi_2-\phi_1}{2}\right)} + \cos{\phi_1}\cos{\phi_2} \sin^2{\left(\frac{\lambda_2-\lambda_1}{2}\right)}, \end{equation}\]

which looks a little simpler. We will work with these forms throughout.

The great circle distance, \(d\), is the shorter arc joining two points on a great circle. We can also consider the chord (straight line) joining the two points, and we let its length be \(C\).

We can immediately observe some relationships between \(d\), \(C\) and the angle \(\sigma\) (measured in radians) that the great circle arc makes with the centre of the sphere: we have \(d=\sigma\) (as the radius is \(1\)), and by drawing the isosceles triangle \(POQ\), we see that \(C=2\sin\bigr(\frac{\sigma}{2}\bigr)\):

Putting these together by eliminating \(\sigma\), \[\frac{C}{2}=\sin\left(\frac{d}{2}\right).\] If we square this and multiply by \(4\), we get an alternative expression: \[\begin{equation} \label{eq:chordhaversine} C^2=4\sin^2\left(\frac{d}{2}\right)=4\mathop{\mathrm{haversin}} d. \end{equation}\]We now offer four approaches to deducing the haversine formula. The first two of these use coordinates, so we begin with this.

#### Coordinates

Given a point on the unit sphere with latitude \(\phi\) and longitude \(\lambda\), we can calculate the corresponding Cartesian coordinates, as shown in the following diagram.

What are the Cartesian coordinates of \(P\) (that is, its \((x,y,z)\) coordinates)?

Once we know the Cartesian coordinates of the points \(P\) and \(Q\) with latitudes \(\phi_1\) and \(\phi_2\) and longitudes \(\lambda_1\) and \(\lambda_2\) respectively, we can calculate the differences in the \(x\), \(y\) and \(z\) coordinates of the two points and use those to calculate appropriate distances or angles in order to find \(d\).

However, the expressions are going to become quite cumbersome with so many variables, so we look for a way to simplify things before going further.

One trick is to notice that the distances and angle \(\sigma\) all remain unchanged if we rotate the sphere about the north-south axis. This only changes the longitudes of the points, so we could rotate it until \(P\) has longitude \(0\) and \(Q\) has longitude \(\lambda_2-\lambda_1\), which we will call \(\Delta\lambda\) (the difference between the \(\lambda\)s).

Then the coordinates of \(P\) and \(Q\) become: \[\begin{align*} P &= (\cos\phi_1, 0, \sin\phi_1)\\ Q &= (\cos\phi_2\cos(\Delta\lambda), \cos\phi_2\sin(\Delta\lambda), \sin\phi_2). \end{align*}\]We could perhaps try rotating further so that \(P\) lies at the north pole, say, but it is not clear what would happen to \(Q\) if we did that, so we will stop rotating at this point.

#### Approach 1: Pythagoras and trigonometry

What is the distance between \(P\) and \(Q\)?

Recalling from above that \(\frac{C}{2}=\sin\bigl(\frac{d}{2}\bigr)\), can you use this to obtain an expression for \(\sin^2\bigl(\frac{d}{2}\bigr)\) as it appears in the haversine formula (equation \(\eqref{eq:haversinesine}\) above)?

#### Approach 2: Vectors

Since we know the coordinates of \(P\) and \(Q\), we also know the vectors \(\mathbf{p}=\overrightarrow{OP}\) and \(\mathbf{q}=\overrightarrow{OQ}\). They each have unit length, so the angle between them, \(\sigma\), is given by \(\cos\sigma=\mathbf{p}\cdot\mathbf{q}\).

Can you use this to find an expression for \(\sin^2\bigl(\frac{d}{2}\bigr)\)?

We can then continue exactly as in Approach 1.

#### Approach 3: Spherical trigonometry

Just as there are sine and cosine rules for triangles in the plane, there are similar rules for triangles on a sphere.

Here is a drawing of a spherical triangle:

The angles \(A\), \(B\) and \(C\) are the angles at the vertices \(A\), \(B\) and \(C\) of the triangle, while \(a\), \(b\) and \(c\) are the side lengths as usual. But if we take the sphere to have radius \(1\), then they can also be thought of as angles: \(a\) is the angle \(\angle BOC\), for example. We therefore assume the sphere to have radius \(1\) from now on.

There are two spherical cosine rules and one spherical sine rule, which are:

spherical cosine rule (1): \[\cos a=\cos b\cos c+\sin b\sin c\cos A\] (and similar equations with \(a\), \(b\), \(c\) and \(A\), \(B\), \(C\) permuted)

spherical cosine rule (2): \[\cos A=-\cos B\cos C+\sin B\sin C\cos a\] (and similar equations with \(a\), \(b\), \(c\) and \(A\), \(B\), \(C\) permuted)

spherical sine rule: \[\dfrac{\sin A}{\sin a}=\dfrac{\sin B}{\sin b}=\dfrac{\sin C}{\sin c}\]

Can you use this information to find an expression for \(\cos\sigma\) and then \(\sin^2\bigl(\frac{d}{2}\bigr)\) in our case?

Again, we can then continue exactly as in Approach 1.

#### Approach 4: Using haversines directly

The above approaches have all proven the haversine formula via the familiar sine and cosine functions. In this approach, we will use equation \(\eqref{eq:chordhaversine}\) above for the length of a chord in terms of the haversine function: \[C^2=4\mathop{\mathrm{haversin}} d,\] where \(d\) is the great circle distance (or angle subtended) and \(C\) is the corresponding chord length.

We wish to find the great-circle distance \(d\) between the points \(P\) and \(Q\) with latitudes \(\phi_1\) and \(\phi_2\) and longitudes \(\lambda_1\) and \(\lambda_2\) respectively.

By thinking about symmetry, we are prompted to also consider the points \(R\) and \(S\) with latitudes \(\phi_1\) and \(\phi_2\) but with swapped longitudes \(\lambda_2\) and \(\lambda_1\) respectively. This gives us the following configuration:

In this diagram, the blue lines are all straight chords (though an optical illusion makes it look otherwise), and the red line is the great circle between \(P\) and \(Q\).

Can you use this diagram to deduce the haversine formula?

How do these four approaches compare?

#### Why use haversines?

Surely we could just use \(\cos\bigl(\frac{d}{R}\bigr)=\cos\phi_1\cos\phi_2\cos(\Delta\lambda)+\sin\phi_1\sin\phi_2\) to calculate great circle distances, as we found in Approaches 2 and 3? Why would anyone need the haversine? (We have reinserted the radius \(R\) here, as practical calculations will use standard units of measurement. Also, a reason is given in Lost but lovely: the haversine for introducing the haversine, but that was just to simplify the formula involving three occurrences of \(\sin^2\frac{x}{2}\).)

Navigators in the past would have had to make do with tables and slide rules: there were no calculators. So what would happen if we try to work out the distance between Dover and Calais, say, using the two different formulae, but restricting ourselves to the limited calculating tools of the past?

Here is some relevant information:

- Dover: latitude \(\quantity{51.15^\circ}{N}\), longitude \(\quantity{1.33^\circ}{E}\)
- Calais: latitude \(\quantity{50.97^\circ}{N}\), longitude \(\quantity{1.85^\circ}{E}\)
- Radius of the Earth: \(\quantity{6378}{km}\)

And trigonometric values, accurate to 4 significant figures or 4 decimal places:

\[\begin{align*} \sin 51.15^\circ&=0.7788 &\cos 51.15^\circ&=0.6273 \\ \sin 50.97^\circ&=0.7768 &\cos 50.97^\circ&=0.6297 \\ \cos 0.52^\circ&=1.0000\\ \mathop{\mathrm{haversin}} 0.18^\circ&=2.467\times10^{-6}\\ \mathop{\mathrm{haversin}} 0.52^\circ&=2.059\times10^{-5} \end{align*}\]Using these values, calculate the great circle distance from Calais to Dover using the two different formulae (and you can use a calculator for the multiplications, and for working out the final inverse cosine and haversine). To truly emulate the past, round all of your intermediate calculations to 4 significant figures as well.