This page presents a variety of calculations for latitude/longitude points, with the formulæ and code fragments for implementing them.
All these formulæ are for calculations on the basis of a spherical earth (ignoring ellipsoidal effects) – which is accurate enough^{*} for most purposes… [In fact, the earth is very slightly ellipsoidal; using a spherical model gives errors typically up to 0.3% – see notes for further details].
And you can see it on a map (aren’t those Google guys wonderful!)
This uses the ‘haversine’ formula to calculate the greatcircle distance between two points – that is, the shortest distance over the earth’s surface – giving an ‘asthecrowflies’ distance between the points (ignoring any hills, of course!).
Haversine formula: 
a = sin²(Δlat/2) + cos(lat_{1}).cos(lat_{2}).sin²(Δlong/2) c = 2.atan2(√a, √(1−a)) d = R.c 
where R is earth’s radius (mean radius = 6,371km); note that angles need to be in radians to pass to trig functions! 

JavaScript:  var R = 6371; // km var dLat = (lat2lat1).toRad(); var dLon = (lon2lon1).toRad(); var lat1 = lat1.toRad(); var lat2 = lat2.toRad(); var a = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1a)); var d = R * c; 
The haversine formula^{1} ‘remains particularly wellconditioned for numerical computation even at small distances’ – unlike calculations based on the spherical law of cosines. The ‘versed sine’ is 1cosθ, and the ‘halfversedsine’ (1cosθ)/2 = sin²(θ/2) as used above. It was published by R W Sinnott in Sky and Telescope, 1984, though known about for much longer by navigators. (For the curious, c is the angular distance in radians, and a is the square of half the chord length between the points). A (surprisingly marginal) performance improvement can be obtained, of course, by factoring out the terms which get squared.
In fact, when Sinnott published the haversine formula, computational precision was limited. Nowadays, JavaScript (and most modern computers & languages) use IEEE 754 64bit floatingpoint numbers, which provide 15 significant figures of precision. With this precision, the simple spherical law of cosines formula (cos c = cos a cos b + sin a sin b cos C) gives wellconditioned results down to distances as small as around 1 metre. (Note that the geodetic form of the law of cosines is rearranged from the canonical one so that the latitude can be used directly, rather than the colatitude).
This makes the simpler law of cosines a reasonable 1line alternative to the haversine formula for many purposes. The choice may be driven by coding context, available trig functions (in different languages), etc.
Spherical law of cosines: 
d = acos(sin(lat_{1}).sin(lat_{2})+cos(lat_{1}).cos(lat_{2}).cos(long_{2}−long_{1})).R 
JavaScript:  var R = 6371; // km var d = Math.acos(Math.sin(lat1)*Math.sin(lat2) + Math.cos(lat1)*Math.cos(lat2) * Math.cos(lon2lon1)) * R; 
Excel:  =ACOS(SIN(lat1)*SIN(lat2)+COS(lat1)*COS(lat2)*COS(lon2lon1))*6371 
(Note that here and in all subsequent code fragments, for simplicity I do not show conversions from degrees to radians; see below for complete versions).
If performance is an issue and accuracy less important, for small distances Pythagoras’ theorem can be used on an equirectangular projection:^{*}
x = Δlon.cos(lat) y = Δlat d = R.√x² + y² 

JavaScript:  var x = (lon2lon1) * Math.cos((lat1+lat2)/2);
var y = (lat2lat1);
var d = Math.sqrt(x*x + y*y) * R;
(lat/lon in radians!) 
This uses just one trig and one sqrt function – as against halfadozen trig functions for cos law, and 7 trigs + 2 sqrts for haversine. Accuracy is somewhat complex: along meridians there are no errors, otherwise they depend on distance, bearing, and latitude, but are small enough for many purposes^{*} (and often trivial compared with the spherical approximation itself).
In general, your current heading will vary as you follow a great circle path (orthodrome); the final heading will differ from the initial heading by varying degrees according to distance and latitude (if you were to go from say 35°N,45°E (Baghdad) to 35°N,135°E (Osaka), you would start on a heading of 60° and end up on a heading of 120°!).
This formula is for the initial bearing (sometimes referred to as forward azimuth) which if followed in a straight line along a greatcircle arc will take you from the start point to the end point:^{1}
Formula:  θ =  atan2(  sin(Δlong).cos(lat_{2}), cos(lat_{1}).sin(lat_{2}) − sin(lat_{1}).cos(lat_{2}).cos(Δlong) ) 
JavaScript:  var y = Math.sin(dLon) * Math.cos(lat2); var x = Math.cos(lat1)*Math.sin(lat2)  Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon); var brng = Math.atan2(y, x).toDeg(); 

Excel:  =ATAN2(COS(lat1)*SIN(lat2)SIN(lat1)*COS(lat2)*COS(lon2lon1), SIN(lon2lon1)*COS(lat2)) * Note that Excel reverses the arguments to ATAN2 – see notes below 
Since atan2 returns values in the range π ... +π (that is, 180° ... +180°), to normalise the result to a compass bearing (in the range 0° ... 360°, with ve values transformed into the range 180° ... 360°), convert to degrees and then use (θ+360) % 360, where % is modulo.
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
This is the halfway point along a great circle path between the two points.^{1}
Formula:  Bx = cos(lat_{2}).cos(Δlong) By = cos(lat_{2}).sin(Δlong) lat_{m} = atan2(sin(lat_{1}) + sin(lat_{2}), √((cos(lat_{1})+Bx)² + By²)) lon_{m} = lon_{1} + atan2(By, cos(lat_{1})+Bx) 
JavaScript:  var Bx = Math.cos(lat2) * Math.cos(dLon); var By = Math.cos(lat2) * Math.sin(dLon); var lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2), Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) ); var lon3 = lon1 + Math.atan2(By, Math.cos(lat1) + Bx); 
Just as the initial bearing may vary from the final bearing, the midpoint may not be located halfway between latitudes/longitudes; the midpoint between 35°N,45°E and 35°N,135°E is around 45°N,90°E.
Given a start point, initial bearing, and distance, this will calculate the destination point and final bearing travelling along a (shortest distance) great circle arc.
Formula:  lat_{2} = asin(sin(lat_{1})*cos(d/R) + cos(lat_{1})*sin(d/R)*cos(θ)) 
lon_{2} = lon_{1} + atan2(sin(θ)*sin(d/R)*cos(lat_{1}), cos(d/R)−sin(lat_{1})*sin(lat_{2}))  
θ is the bearing (in radians, clockwise from
north); d/R is the angular distance (in radians), where d is the distance travelled and R is the earth’s radius 

JavaScript:  var lat2 = Math.asin( Math.sin(lat1)*Math.cos(d/R) + Math.cos(lat1)*Math.sin(d/R)*Math.cos(brng) ); var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(lat1), Math.cos(d/R)Math.sin(lat1)*Math.sin(lat2)); 
Excel:  lat2: =ASIN(SIN(lat1)*COS(d/R) + COS(lat1)*SIN(d/R)*COS(brng)) lon2: =lon1 + ATAN2(COS(d/R)SIN(lat1)*SIN(lat2), SIN(brng)*SIN(d/R)*COS(lat1)) 
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
This is a rather more complex calculation than most others on this page, but I've been asked for it a number of times. See below for the JavaScript.
Formula:  d_{12} = 2.asin( √(sin²(Δlat/2)
+ cos(lat_{1}).cos(lat_{2}).sin²(Δlon/2)) ) if sin(lon_{2}−lon_{1}) > 0 α_{1} = (θ_{1} − θ_{12} + π) % 2.π − π α_{3} = acos( −cos(α_{1}).cos(α_{2}) + sin(α_{1}).sin(α_{2}).cos(d_{12})
) 
where  lat_{1}, lon_{1}, θ_{1} : 1st point & bearing % = mod 
note –  if sin(α_{1})=0 and sin(α_{2})=0: infinite solutions if sin(α_{1}).sin(α_{2}) < 0: ambiguous solution this formulation is not always wellconditioned for meridional or equatorial lines 
Note this can also be solved using vectors rather than trigonometry:
Here’s a new one: I’ve sometimes been asked about distance of a point from a greatcircle path (sometimes called cross track error).
Formula:  d_{xt} = asin(sin(d_{13}/R)*sin(θ_{13}−θ_{12})) * R 
where  d_{13} is distance from start point to third point θ_{13} is (initial) bearing from start point to third point θ_{12} is (initial) bearing from start point to end point R is the earth’s radius 
JavaScript:  var dXt = Math.asin(Math.sin(d13/R)*Math.sin(brng13brng12)) * R; 
Here, the greatcircle path is identified by a start point and an end point – depending on what initial data you’re working from, you can use the formulæ above to obtain the relevant distance and bearings. The sign of d_{xt} tells you which side of the path the third point is on.
The alongtrack distance, from the start point to the closest point on the path to the third point, is
Formula:  d_{at} = acos(cos(d_{13}/R)/cos(d_{xt}/R)) * R 
where  d_{13} is distance from start point to third point d_{xt} is crosstrack distance R is the earth’s radius 
JavaScript:  var dAt = Math.acos(Math.cos(d13/R)/Math.cos(dXt/R)) * R; 
And: ‘Clairaut’s formula’ will give you the maximum latitude of a great circle path, given a bearing and latitude on the great circle:
Formula:  lat_{max} = acos(abs(sin(θ)*cos(lat))) 
JavaScript:  var latMax = Math.acos(Math.abs(Math.sin(brng)*Math.cos(lat))); 
A ‘rhumb line’ (or loxodrome) is a path of constant bearing, which crosses all meridians at the same angle.
Sailors used to (and sometimes still) navigate along rhumb lines since it is easier to follow a constant compass bearing than to be continually adjusting the bearing, as is needed to follow a great circle. Rhumb lines are straight lines on a Mercator Projection map (also helpful for navigation).
Rhumb lines are generally longer than greatcircle (orthodrome) routes. For instance, London to New York is 4% longer along a rhumb line than along a great circle – important for aviation fuel, but not particularly to sailing vessels. New York to Beijing – close to the most extreme example possible (though not sailable!) – is 30% longer along a rhumb line.
No, I’ve not included decimal minutes: a decimal system is easy, a sexagesimal system has merits, but mixing the two is a complete sow’s ear. Switch off the option on your GPS!
If you have any queries or find any problems, contact me at [email protected].
© 20022011 Chris Veness