Vincenty’s formula as it is used in the script:
| a, b = major & minor semiaxes of the ellipsoid | ||
| f = flattening (a−b)/a | ||
| φ1, φ2 = geodetic latitude | ||
| L = difference in longitude | ||
| U1 = atan((1−f).tanφ1) (U is ‘reduced latitude’) | ||
| U2 = atan((1−f).tanφ2) | ||
| λ = L (first approximation) | ||
| iterate until change in λ is negligible (e.g. 10-12 ≈ 0.06mm) { | ||
| sinσ = √[ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ] | (14) | |
| cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ | (15) | |
| σ = atan2(sinσ, cosσ) | (16) | |
| sinα = cosU1.cosU2.sinλ / sinσ | (17) | |
| cos²α = 1 − sin²α (trig identity; §6) | ||
| cos2σm = cosσ − 2.sinU1.sinU2/cos²α | (18) | |
| C = f/16.cos²α.[4+f.(4−3.cos²α)] | (10) | |
| λ′ = L + (1−C).f.sinα.{σ+C.sinσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]} | (11) | |
| } | ||
| u² = cos²α.(a²−b²)/b² | ||
| A = 1+u²/16384.{4096+u².[−768+u².(320−175.u²)]} | (3) | |
| B = u²/1024.{256+u².[−128+u².(74−47.u²)]} | (4) | |
| Δσ = B.sinσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)]} | (6) | |
| s = b.A.(σ−Δσ) | (19) | |
| α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ) | (20) | |
| α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ) | (21) | |
Where:
- s is the distance (in the same units as a & b)
- α1 is the initial bearing, or forward azimuth
- α2 is the final bearing (in direction p1→p2)
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, and GRS-67 in South America. America (NAD83) and Australia (GDA) use GRS-80, functionally equivalent to the WGS-84 ellipsoid.
| WGS-84 | a = 6 378 137 m (±2 m) | b = 6 356 752.3142 m | f = 1 / 298.257223563 | |
| GRS-80 | a = 6 378 137 m | b = 6 356 752.3141 m | f = 1 / 298.257222101 | |
| Airy (1830) | a = 6 377 563.396 m | b = 6 356 256.909 m | f = 1 / 299.3249646 | |
| Int’l 1924 | a = 6 378 388 m | b = 6 356 911.946 m | f = 1 / 297 | |
| Clarke (1880) | a = 6 378 249.145 m | b = 6 356 514.86955 m | f = 1 / 293.465 | |
| GRS-67 | a = 6 378 160 m | b = 6 356 774.719 m | f = 1 / 298.25 |
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.
Some of the terms involved are explained in Ed Williams’ notes on Spheroid Geometry.
| Test case (from Geoscience Australia), using WGS-84: | ||
| Flinders Peak | 37°57′03.72030″S, 144°25′29.52440″E | |
| Buninyong | 37°39′10.15610″S, 143°55′35.38390″E | |
| s | 54 972.271 m | |
| α1 | 306°52′05.37″ | |
| α2 | 127°10′25.07″ (≡ 307°10′25.07″ p1→p2) | |
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 LGPL 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, 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-2009 Chris Veness
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2009 */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
* Calculate geodesic distance (in m) between two points specified by latitude/longitude
* (in numeric degrees) using Vincenty inverse formula for ellipsoids
*/
function distVincenty(lat1, lon1, lat2, lon2) {
var a = 6378137, b = 6356752.3142, f = 1/298.257223563; // WGS-84 ellipsiod
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;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
* extend String object with method for parsing degrees or lat/long values to numeric degrees
*
* this is very flexible on formats, allowing signed decimal degrees, or deg-min-sec suffixed by
* compass direction (NSEW). A variety of separators are accepted (eg 3º 37' 09"W) or fixed-width
* format without separators (eg 0033709W). Seconds and minutes may be omitted. (Minimal validation
* is done).
*/
String.prototype.parseDeg = function() {
if (!isNaN(this)) return Number(this); // signed decimal degrees without NSEW
var degLL = this.replace(/^-/,'').replace(/[NSEW]/i,''); // strip off any sign or compass dir'n
var dms = degLL.split(/[^0-9.]+/); // split out separate d/m/s
for (var i in dms) if (dms[i]=='') dms.splice(i,1); // remove empty elements (see note below)
switch (dms.length) { // convert to decimal degrees...
case 3: // interpret 3-part result as d/m/s
var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600; break;
case 2: // interpret 2-part result as d/m
var deg = dms[0]/1 + dms[1]/60; break;
case 1: // decimal or non-separated dddmmss
if (/[NS]/i.test(this)) degLL = '0' + degLL; // - normalise N/S to 3-digit degrees
var deg = dms[0].slice(0,3)/1 + dms[0].slice(3,5)/60 + dms[0].slice(5)/3600; break;
default: return NaN;
}
if (/^-/.test(this) || /[WS]/i.test(this)) deg = -deg; // take '-', west and south as -ve
return deg;
}
// note: whitespace at start/end will split() into empty elements (except in IE)
/*
* extend Number object with methods for converting degrees/radians
*/
Number.prototype.toRad = function() { // convert degrees to radians
return this * Math.PI / 180;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
<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(lat1.value.parseDeg(), long1.value.parseDeg(),
lat2.value.parseDeg(), long2.value.parseDeg()) + ' m'">
<input name="result" value="">