Great Circle on Spherical Earth

The shortest path between two points on the (spherical) Earth lies on the great circle going through these two points. The great circle is - as the name suggests - the largest possible circle on the sphere. All great circles of a given sphere have both the same radius and the same center as the sphere. The great circle's plane passes through the center of the sphere and divides it into 2 equal hemispheres.

Great Circle between Two Points

The equation of the great circle is computed most easily in cartesian coordinates. So, if the points were given in spherical coordinates we would transform them to cartesian coordinates first.

Let

be the vectors corresponding to the two points P₁ and P₂, respectively. The great circle is located in the plane OP₁P₂ spanned by these two vectors.

If we find two orthonormal vectors u and v in this plane then the equation of the great circle will be

The obvious choice for the vector u is the normalized vector p₁:

The vector product

is orthogonal to the plane OP₁P₂. The unit vector

is in the plane OP₁P₂ and is orthogonal to u.

Another way to compute the vector v is by using the Gram-Schmidt orthonormalization process. The angle δ between p₁ and p₂ can be computed using the scalar product:

The projection of p₂ to u is:

The vector

is orthogonal to p₁ and is in the plane OP₁P₂. Hence the normalized vector v is our searched vector:

The above calculations break if the vectors p₁ and p₂ are parallel, because then they do not specify a plane. This happens if the corresponding points are either identical or antipodal. In this case there are infinite many great circles through P₁ and P₂. Their equations are still given by

where v is an arbitrary unit vector orthogonal to u.

Shortest Distance between Two Points

The shortest path between two points on the earth surface is the circular arc on the great circle between these points. So, the shortest distance d is

where R is the earth radius and δ is the angle between the two points. From the equation (6) it follows:

By substituting spherical coordinates and using the identity:

we get:

Finally, the shortest distance d is:

Heading (Bearing, Azimuth)

Heading (also called bearing or azimuth) is the angle α between a great circle and a meridian. Heading plays an important role in navigation. When we travel from the point P₁ to the point P₂ on the shortest path, then the heading at point P₁ specifies the starting direction of our journey. The heading changes during the voyage, of course.

Heading as angle between planes

Heading α is the angle between the great circle plane OP₁P₂ and the meridian plane OP₁P₂, or, equivalently between the respective normal vectors of these two planes. The vector orthogonal to the meridian plane is

The vector orthogonal to the path plane is w = (w_x, w_y, w_z)^T = p₁ x p₁.

The dot product of the normalized vectors w and m equals the cosine of the heading angle α:

When we can take arccos() we get the angle α, 0 ≤ α ≤ 180°. The usual convention is that a positive angle denotes travelling eastwards while a negative angle denotes travelling westwards. Under the assumption that we travel along the shorter arc from P₁ to P₂ we can determine the direction by considering the vector product w. When it points upwards with respect to the positive z-axis the journey goes eastwards, when it points downwards the journey goes westwards. So, our formula for the heading angle is finally:

Simplification of the geometry

To simplify the following calculations I will make two assumptions without loss of generality:

1. The radius of the earth is 1. This is justified because only angles are considered.
2. The point P₁ has the longitude 0 and the point P₂ has the longitude d λ = λ₂ - λ₁. I just start counting the longitude at the meridian passing through P₁. The vector p₁ is hence (x₁, 0, z₁)^T.

Under the simplifying assumptions we get:

and

and

Heading as angle in a spherical triangle

The heading is the angle at the vertex P₁ of the spherical triangle NP₁P₂, where N is the north pole. The sides of this triangle are arcs of the respective great circles and their lengths are equal to the respective angles at the center of the Earth.

Now we can apply the cosine rule formula for the side P₂N that is opposite to the angle α (see the equation (10) in Weisstein, Eric W. Spherical Trigonometry in MathWorld):

Solving for cosα and taking into account that

we get:

Tangens formula

Let us have a look at the geometry of the vector w:

In the equation (21) we computed in fact α using the cosine formula in the triangle OBC. Now, we will compute it using the tangens formula:

Here, w_xz is the projection of the vector w onto the XZ-plane:

The last equality follows from the fact that |p₁| = 1 and y₁ = 0. Substituting (26) and (20b) into (25) we get:

Finally, using spherical coordinates we obtain:

Even if this formula looks more complicated than the formula (24) it is often recommended because it avoids the calculation of the angle δ between the two points P₁ and P₂. However, significant care must be taken by applying atan() to the equation (28). A simpler way is to use atan2() instead of atan():

dl = lon2 - lon1
y2 = sin(dl) * cos(lat2)
wy = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dl)
alpha = atan2(y2, wy)

References

• Wikipedia: Great-circle distance
• Chris Veness: Calculate distance, bearing and more between Latitude/Longitude points
• Weisstein, Eric W. "Great Circle." From MathWorld--A Wolfram Web Resource
• Weisstein, Eric W. "Spherical Trigonometry" From MathWorld--A Wolfram Web Resource
• PROJ.4 - Geodesic Calculations
• Ed Williams: Aviation Formulary
• James R. Clynch: Paths Between Points on Earth
• Mathforum: Latitude & Longitude
• Mathforum: Bearing Between Two Points