Movable Type Scripts

---

Vincenty formula for distance between two Latitude/Longitude points

Vincenty’s formula as it is used in the script:

Where:

• s is the distance (in the same units as a & b)
• α₁ is the initial bearing, or forward azimuth
• α₂ is the final bearing (in direction p₁→p₂)

Note: Vincenty observes that eqn. (18) becomes indeterminate over equatorial lines (since cos²α → 0); in this case, set cos2σ_m to 0 and the result is computed correctly. He also points out that the formula may have no solution between two nearly antipodal points; an iteration limit traps this case (Vincenty says “this will occur when λ, as computed by eqn. (11), is greater than π in absolute value”, but this is not always a reliable test).
Note: some implementations of Vincenty’s formula inefficiently use a large number of trig functions; Vincenty devised this solution with an eye for efficiency in implementation, and this one uses just one each of sin, cos, sqrt, and atan2 for each iteration – only 3 or 4 iterations are generally required. [Formulation updated Dec 05 to make it closer to Vincenty’s original and computationally more efficient.]

The most accurate and widely used globally-applicable model for the earth ellipsoid is WGS-84, used in this script. Other ellipsoids offering a better fit to the local geoid include Airy (1830) in the UK, International 1924 in much of Europe, Clarke (1880) in Africa, GRS-67 in South America, and many others. America (NAD83) and Australia (GDA) use GRS-80, functionally equivalent to the WGS-84 ellipsoid.

Note that to locate latitude/longitude points on these ellipses, they are associated with specific datums: for instance, OSGB36 for Airy in the UK, ED50 for Int’l 1924 in Europe; WGS-84 defines a datum as well as an ellipse. See my convert coordinates page for converting points between different datums.

Aside for fellow geeks: it is often said that WGS-84 and GRS-80 are ‘functionally equivalent’. Well, I would concur that a 0.1 mm difference in Earth radius is pretty equivalent, but why any difference? Apparently, when defining WGS-84, the normalized second degree zonal harmonic gravitational coefficient C_2,0, derived from the GRS-80 value for the dynamic form factor J₂, was truncated to 8 significant digits in the normalization process!

Some of the terms involved are explained in Ed Williams’ notes on Spheroid Geometry.

Notes:

• Trig functions take arguments in radians, so latitude, longitude, and bearings in degrees (either decimal or degrees/minutes/seconds) need to be converted to radians, rad = π.deg/180. When converting radians back to degrees (deg = 180.rad/π), West is negative if using signed decimal degrees. For bearings, values in the range -π to +π [-180° to +180°] need to be converted to 0 to +2π [0°–360°]; this can be done by (brng+2.π)%2.π [brng+360)%360] where % is the modulo operator.
• The atan2() function used here takes two arguments, atan2(y, x), and computes the arc tangent of the ratio y/x. It is more flexible than atan(y/x), since it handles x=0, and it also returns values in all 4 quadrants -π to +π (the atan function returns values in the range -π/2 to +π/2).
• If you implement any formula involving atan2 in Microsoft Excel, you will need to reverse the arguments, as Excel has them the opposite way around from JavaScript – conventional order is atan2(y, x), but Excel uses atan2(x, y). To use atan2 in a (VBA) macro, you can use WorksheetFunction.Atan2(). You will also have to rename the variables A, B / a, b, as VBA is case-insensitive.
• All bearings are with respect to true north, 0°=N, 90°=E, etc; if you are working from a compass, magnetic north varies from true north in a complex way around the earth, and the difference has to be compensated for by variances indicated on local maps.
• For miles, divide km by 1.609344
• For nautical miles, divide km by 1.852

See below for the source code of the JavaScript implementation. These functions should be simple to translate into other languages if required.

I offer these formulæ & scripts for free use and adaptation as my contribution to the open-source info-sphere from which I have received so much. You are welcome to re-use these scripts [under a simple attribution license, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.

If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.

© 2002-2011 Chris Veness

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2011             */
/*                                                                                                */
/* from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
/*       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975    */
/*       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf                                             */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

/**
 * Calculates geodetic distance between two points specified by latitude/longitude using 
 * Vincenty inverse formula for ellipsoids
 *
 * @param   {Number} lat1, lon1: first point in decimal degrees
 * @param   {Number} lat2, lon2: second point in decimal degrees
 * @returns (Number} distance in metres between points
 */
function distVincenty(lat1, lon1, lat2, lon2) {
  var a = 6378137, b = 6356752.314245,  f = 1/298.257223563;  // WGS-84 ellipsoid params
  var L = (lon2-lon1).toRad();
  var U1 = Math.atan((1-f) * Math.tan(lat1.toRad()));
  var U2 = Math.atan((1-f) * Math.tan(lat2.toRad()));
  var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
  var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
  
  var lambda = L, lambdaP, iterLimit = 100;
  do {
    var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
    var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
      (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
    if (sinSigma==0) return 0;  // co-incident points
    var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
    var sigma = Math.atan2(sinSigma, cosSigma);
    var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    var cosSqAlpha = 1 - sinAlpha*sinAlpha;
    var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    if (isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
    var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
    lambdaP = lambda;
    lambda = L + (1-C) * f * sinAlpha *
      (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
  } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);

  if (iterLimit==0) return NaN  // formula failed to converge

  var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
  var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
    B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
  var s = b*A*(sigma-deltaSigma);
  
  s = s.toFixed(3); // round to 1mm precision
  return s;
  
  // note: to return initial/final bearings in addition to distance, use something like:
  var fwdAz = Math.atan2(cosU2*sinLambda,  cosU1*sinU2-sinU1*cosU2*cosLambda);
  var revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
  return { distance: s, initialBearing: fwdAz.toDeg(), finalBearing: revAz.toDeg() };
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
  

Usage:

<p>Lat 1: <input name="lat1" value="53 09 02N"> Long 1: <input name="long1" value="001 50 40W"></p>
<p>Lat 2: <input name="lat2" value="52 12 19N"> Long 2: <input name="long2" value="000 08 33W"></p>

<input type="button" value="calculate distance" 
    onClick="result.value = distVincenty(Geo.parseDMS(lat1.value), Geo.parseDMS(long1.value), 
                                         Geo.parseDMS(lat2.value), Geo.parseDMS(long2.value)) + ' m'">
<input name="result" value="">
  

See the Lat/Long page for the deg/min/sec parsing method Geo.parseDMS()