### By Ed Williams

#### Introduction

This introduction is written for pilots (and others) who are interested in great circle navigation and would like to know how to compute courses, headings and other quantities of interest. These formulae can be programmed into your calculator or spreadsheet. I’ll attempt to include enough information that those familiar with plane trigonometry can derive additional results if required.

It is a well known that the shortest distance between two points is a straight line. However anyone attempting to fly from Los Angeles to New York on the straight line connecting them would have to dig a very substantial tunnel first. The shortest distance, *following the earth’s surface *lies vertically above the aforementioned straight line route. This route can be constructed by slicing the earth in half with an imaginary plane through LAX and JFK. This plane cuts the (assumed spherical) earth in a circular arc connecting the two points, called a *great *circle. Only planes through the center of the earth give rise to great circles. Any plane will cut a sphere in a circle, but the resulting little circles are not the shortest distance between the points they connect. A little thought will show that lines of longitude (meridians) are great circles, but lines of latitude, with the exception of the equator, are not.

I will assume the reader is familiar with latitude and longitude as a means of designating locations on the earth’s surface. For the convenience of North Americans I will take North latitudes and West longitudes as positive and South and East negative. The longitude is the opposite of the usual mathematical convention. True course is defined as usual, as the angle between the course line and the local meridian measured clockwise.

The first important fact to realise is that in general a great circle route has a true course that varies from point to point. For instance the great circle route between two points of equal (non-zero) latitude does not follow the line of latitude in an E-W direction, but arcs towards the pole. It *is* possible to fly between two points using an unvarying true course, but in general the resulting route differs from the great circle route and is called a *rhumb *line. Unlike a great circle which encircles the earth in a closed curve, a pilot flying a rhumb line would spiral poleward.

Natural questions are to seek the great circle distance between two specified points and true course at points along the route. The required spherical trigonometric formulae are greatly simplified if angles and distances are measured in the appropriate natural units, which are both radians! A radian, by definition, is the angle subtended by a circular arc of unit length and unit radius. Since the length of a complete circular arc of unit radius is 2*pi, the conversion is 360 degrees equals 2*pi radians, or:

angle_radians=(pi/180)*angle_degrees angle_degrees=(180/pi)*angle_radians

Great circle distance can be likewise be expressed in radians by defining the distance to be the angle subtended by the arc at the center of the earth. Since by definition, one nautical mile subtends one minute (=1/60 degree) of arc, we have:

distance_radians=(pi/(180*60))*distance_nm distance_nm=((180*60)/pi)*distance_radians

In all subsequent formulae all distances and angles, such as latitudes, longitudes and true courses will be assumed to be given in radians, greatly simplifying them, and in applications the above formulae and their inverses are necessary to convert back and forth between natural and practical units. Examples of this process are given later.

Note: the nautical mile is currently defined to be 1852 meters – which to be consistent with its historical definition implies the earth’s radius to be 1.852 * (180*60/pi) = 6366.71 km, which indeed lies between the currently accepted ( WGS84) equatorial and polar radii of 6378.137 and 6356.752 km, respectively. Other choices of the earth’s radius in this range are consistent with the spherical approximation and may for some specialized purposes be preferred.

For instance, the former FAI standard for aviation records used the “FAI sphere” with radius 6371.0 km.

To use a different spherical radius, use

distance_km = radius_km * distance_radians distance_radians = distance_km/radius_km

to convert between angular and distance units.

A sample implementation of many of these formulae in the form of an Excel spreadsheet can be found here. The formulae are in VBA macros, for readability, so you will need to enable macros to have them work. If you are unable to open the spreadsheet with macros disabled (to check for viruses) etc, then you may need to patch your Excel. Try http://office.microsoft.com/downloaddetails/Xl8p9pkg.htm for Excel 97 SR-2.

*If you decide to program up these formulae, you’d be well-advised to look at the implementation notes and check your results against the worked examples and spreadsheets.*

#### Some great circle formulae:

##### Distance between points

The great circle distance d between two points with coordinates {lat1,lon1} and {lat2,lon2} is given by:

d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))

A mathematically equivalent formula, which is less subject to rounding error for short distances is:

d=2*asin(sqrt((sin((lat1-lat2)/2))^2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))

##### Lat/lon given radial and distance

A point {lat,lon} is a distance d out on the tc radial from point 1 if:

lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc)) IF (cos(lat)=0) lon=lon1 // endpoint a pole ELSE lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi ENDIF

This algorithm is limited to distances such that dlon <pi/2, i.e those that extend around less than one quarter of the circumference of the earth in longitude. A completely general, but more complicated algorithm is necessary if greater distances are allowed:

lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc)) dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat)) lon=mod( lon1-dlon +pi,2*pi )-pi

#### Implementation notes:

##### Notes on mathematical functions

Note: ^ denotes the exponentiation operator, sqrt is the square root function, acos the arc-cosine (or inverse cosine) function and asin is the arc-sine function. If asin or acos are unavailable they can be implemented using the atan2 function:

acos(x)=atan2(sqrt(1-x^2),x) acos returns a value in the range 0 <= acos <= pi asin(x)=atan2(x,sqrt(1-x^2))} asin returns a value in the range -pi/2 <= asin <= pi/2

Note: Here atan2 has the conventional (C) ordering of arguments, namely atan2(y,x). This is not universal, Excel for instance uses atan2(x,y), but it has asin and acos anyway. Be warned. It returns a value in the range -pi < atan2 <= pi.

Further note: if your calculator/programming language is so impoverished that only atan is available then use:

asin(x)=2*atan(x/(1+sqrt(1-x*x))) acos(x)=2*atan(sqrt((1-x)/(1+x))) x>=0 =pi - 2*atan(sqrt((1+x)/(1-x))) x<0 atan2(y,x)=atan(y/x) x>0 atan2(y,x)=atan(y/x)+pi x<0, y>=0 atan2(y,x)=pi/2 x=0, y>0 atan2(y,x)=atan(y/x)-pi x<0, y<0 atan2(y,x)=-pi/2 x=0, y<0 atan2(0,0) is undefined and should give an error.

Another potential implementation problem is that the arguments of asin and/or acos may, because of rounding error, exceed one in magnitude. With perfect arithmetic this can’t happen. You may need to use “safe” versions of asin and acos on the lines of:

asin_safe(x)=asin(max(-1,min(x,1))) acos_safe(x)=acos(max(-1,min(x,1)))

Note on the mod function. This appears to be implemented differently in different languages, with differing conventions on whether the sign of the result follows the sign of the divisor or the dividend. (We want the sign to follow the divisor or be Euclidean. C’s fmod and Java’s % do not work.) In this document, Mod(y,x) is the remainder on dividing y by x and always lies in the range 0 <=mod <x. For instance: mod(2.3,2.)=0.3 and mod(-2.3,2.)=1.7

If you have a floor function (int in Excel), that returns floor(x)= “largest integer less than or equal to x” e.g. floor(-2.3)=-3 and floor(2.3) =2

mod(y,x) = y - x*floor(y/x) The following should work in the absence of a floor function- regardless of whether "int" truncates or rounds downward: mod=y - x * int(y/x) if ( mod < 0) mod = mod + x

##### Sign Convention

As stated in the introduction, North latitudes and *West *longitudes are treated as positive, and South latitudes and East longitudes negative. It’s easier to go with the flow, but if you prefer another convention you can change the signs in the formulae.

**Reference**

http://williams.best.vwh.net/avform.htm#Intro