# 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

**p**= (x

_{1}_{1}, y

_{1}, z

_{1})

^{T}

**p**= (x

_{1}_{1}, y

_{1}, z

_{1})

^{T}

**p**= (x

_{2}_{2}, y

_{2}, z

_{2})

^{T}

be the vectors corresponding to the two points
*P _{1}* and

*P*, respectively. The great circle is located in the plane

_{2}*OP*spanned by these two vectors.

_{1}P_{2}If we find two orthonormal vectors * u* and

*in this plane then the equation of the great circle will be*

**v****u**cos ω +

**v**sin ω)

The obvious choice for the vector * u* is the normalized vector

*:*

**p**_{1}**u**=

**p**/ R .

_{1}The vector product

**w**=

**p**x

_{1}**p**

_{2}is orthogonal to the plane *OP _{1}P_{2}*.
The unit vector

**v**=

**u**x

**w**/ |

**w**|

is in the plane *OP _{1}P_{2}* and is orthogonal to

*.*

**u**Another way to compute the vector * v* is by using the
Gram-Schmidt orthonormalization process. The angle

*δ*between

*and*

**p**_{1}*can be computed using the scalar product:*

**p**_{2}**p**·

_{1}**p**= R

_{2}^{2}cos δ .

The projection of * p_{2}* to

*is:*

**u****p**= R cos(δ)

_{u}**u**.

The vector

**p**=

_{v}**p**-

_{2}**p**

_{u}is orthogonal to * p_{1}* and is in the plane

*OP*. Hence the normalized vector

_{1}P_{2}*is our searched vector:*

**v****v**=

**p**/ |

_{v}**p**|

_{v}The above calculations break if the vectors * p_{1}* and

*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**_{2}*P*and

_{1}*P*. Their equations are still given by

_{2}**c**= r (

**u**cos ω +

**v**sin ω)

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:

**p**·

_{1}**p**/ R

_{2}^{2}= (x

_{1}x

_{2}+ y

_{1}y

_{2}+ z

_{1}z

_{2}) / R

^{2}.

By substituting spherical coordinates and using the identity:

_{1}- λ

_{2)}= cosλ

_{1}cosλ

_{2}+ sinλ

_{1}sinλ

_{2}

we get:

_{1}- λ

_{2)}cos φ

_{1}cos φ

_{2}+ sin φ

_{1}sin φ

_{2}

Finally, the shortest distance *d* is:

_{1}- λ

_{2)}cos φ

_{1}cos φ

_{2}+ sin φ

_{1}sin φ

_{2}]

## 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 _{1}* to the point

*P*on the shortest path, then the heading at point

_{2}*P*specifies the starting direction of our journey. The heading changes during the voyage, of course.

_{1}### Heading as angle between planes

Heading *α* is the angle between
the great circle plane *OP _{1}P_{2}* and
the meridian plane

*OP*, or, equivalently between the respective normal vectors of these two planes. The vector orthogonal to the meridian plane is

_{1}P_{2}**m**=

**p**x (0, 0, 1)

_{1}^{T}= (y

_{1}, -x

_{1}, 0)

^{T}

The vector orthogonal to the path plane is
* w* =

*(w*=

_{x}, w_{y}, w_{z})^{T}*.*

**p**_{1}x p_{1}The dot product of the normalized vectors * w* and

*equals the cosine of the heading angle*

**m***α*:

**w**/ |

**w**| ·

**m**/ |

**m**|

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 _{1}* to

*P*we can determine the direction by considering the vector product

_{2}*. 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:*

**w**_{z}≥ 0

_{z}< 0

### Simplification of the geometry

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

- The radius of the earth is 1. This is justified because only angles are considered.
- The point
*P*has the longitude 0 and the point_{1}*P*has the longitude_{2}*d λ*=*λ*-_{2}*λ*. I just start counting the longitude at the meridian passing through_{1}*P*. The vector_{1}is hence**p**_{1}*(x*._{1}, 0, z_{1})^{T}

Under the simplifying assumptions we get:

**m**= (0, -x

_{1}, 0)

^{T}

and

_{x}= z

_{1}y

_{2}

_{y}= -x

_{1}z

_{2}+ z

_{1}x

_{2}

_{z}= x

_{1}y

_{2}

and

_{y}/ |

**w**|

### Heading as angle in a spherical triangle

The heading is the angle at the vertex *P _{1}* of the
spherical triangle

*NP*, where

_{1}P_{2}*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 _{2}N*
that is opposite to the angle

*α*(see the equation (10) in Weisstein, Eric W. Spherical Trigonometry in MathWorld):

_{2}= cos Φ

_{1}cos δ + sin Φ

_{1}sin δ cos α

Solving for cos*α* and taking into account that

_{1}= 90° - φ

_{1}

_{2}= 90° - φ

_{2}

we get:

_{2}- sin φ

_{1}cos δ

cos φ

_{1}sin δ

### 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:

_{xz}/ w

_{y}

Here, *w _{xz}* is the projection of the vector

*onto the XZ-plane:*

**w**_{xz}

^{2}= w

_{x}

^{2}+ w

_{z}

^{2}= z

_{1}

^{2}y

_{2}

^{2}+ x

_{1}

^{2}y

_{2}

^{2}

_{xz}

^{2}= y

_{2}

^{2}(x

_{1}

^{2}+ z

_{1}

^{2}) = y

_{2}

^{2}

The last equality follows from the fact that * |p_{1}|* = 1
and

*y*= 0. Substituting (26) and (20b) into (25) we get:

_{1}_{2}

x

_{1}z

_{2}- z

_{1}x

_{2}

Finally, using spherical coordinates we obtain:

_{2}sin dλ

cos φ

_{1}sin φ

_{2}- sin φ

_{1}cos φ

_{2}cos dλ

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 _{1}* 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():

_{2}```
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