 # Movable Type Scripts

## Vincenty solutions of geodesics on the ellipsoid

### Distances and bearings between points on an ellipsoidal-model earth

For the benefit of the terminally obsessive (as well as the genuinely needy, of course), Thaddeus Vincenty devised formulæ for calculating geodesic distances between a pair of lati­tude/longi­tude points on the earth’s surface, using an accurate ellipsoidal model of the earth.

When I looked up the ref­erences (‘Direct and Inverse Solutions of Geodesics on the Ellipsoid with applica­tion of nested equa­tions’), I discovered to my surprise that while the mathematics is utterly beyond me, it is actually quite simple to program.

Vincenty’s solution for the distance between points on an ellipsoidal earth model is accurate to within 0.5 mm distance (!), 0.000015″ bearing, on the ellipsoid being used. Calcula­tions based on a spherical earth model, such as the (much simpler) Haversine, are accurate to around 0.3% – which is still good enough for many (most?) purposes, of course.

Vincenty’s inverse solution can fail on nearly antipodal points. Testing with Geographic­Lib test data, I’ve found this can happen with distances greater than 19,936 km, or within around 75 km of the antipodal point.

Vincenty wrote this at a time when computing resources were expen­sive, and made his solu­tions very computa­tionally efficient – the distance calcula­tion takes around 5 – 10 micro­seconds (hence around 100,000 – 200,000 per second) using Chrome on an aging Core i5 PC.

### Distance/bearing between two points (inverse solution)

Vincenty’s paper presents a set of formulæ rather than a single program; this shows those formulæ pulled together as they are used in the script:

 a, b = major & minor semi-axes of the ellipsoid f = flattening (a−b)/a φ1, φ2 = geodetic latitude L = difference in longitude tan U1/2 = (1−f) · tan φ1/2 (U is ‘reduced latitude’) cos U1/2 = 1 / √1 + tan² U1/2, sin U1/2 = tan U1/2 · cos U1/2 (trig identities; §6) λ = L (first approximation) iterate until change in λ is negligible (e.g. 10-12 ≈ 0.006mm) { sin σ = √[ (cos U2 · sinλ)² + (cos U1 · sin U2 − sin U1 · cos U2 · cos λ)² ] (14) cos σ = sin U1 · sin U2 + cos U1 · cos U2 · cos λ (15) σ = atan(sin σ / cos σ) (16) sin α = cos U1 · cos U2 · sin λ / sin σ (17) cos² α = 1 − sin² α (trig identity; §6) cos 2σm = cos σ − 2 · sin U1 · sin U2 / cos² α (18) C = f/16 · cos² α · [4 + f · (4 − 3 · cos² α)] (10) λ′ = L + (1 − C) · f · sin α · {σ + C · sin σ · [cos 2σ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 σ · {cos 2σm + B/4 · [cos σ · (−1 + 2 · cos² 2σm)     − B / 6 · cos 2σm · (−3 + 4 · sin² σ) · (−3 + 4 · cos² 2σm)]} (6) s = b·A·(σ−Δσ) (19) α1 = atan(cos U2 · sin λ / cos U1 · sin U2 − sin U1 · cos U2 · cos λ) (20) α2 = atan(cos U1 · sin λ / −sin U1 · cos U2 + cos U1 · sin U2 · cos λ) (21)

Where:

• s is the geodesic distance along the surface of the ellipsoid (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 solu­tion between two nearly antipodal points; an itera­tion limit traps this case (Vincenty says “this will occur when λ, as computed by eqn. (11), is greater than π in absolute value”, but in a small number of cases, the formula will not converge even when λ stays below π).

Note: some implementations of Vincenty’s formula inefficiently use a large number of trig functions; Vincenty devised this solu­tion with an eye for efficiency in implementa­tion, and this one uses just one each of sin, cos, sqrt, and atan2 for each itera­tion – only 3 or 4 itera­tions are generally required. Pythagorean trigonometric identities sin²φ + cos²φ = 1 and 1 + tan²φ = sec²φ are used to reduce the number of trig functions invoked.

#### JavaScript

(slightly simplified, ignoring convergence, equatorial lines, coincident/antipodal cases etc)
```const L = λ2 - λ1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanφ.
const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
const tanU2 = (1-f) * Math.tan(φ2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2;

let λ = L, sinλ = null, cosλ = null;    // λ = difference in longitude on an auxiliary sphere
let σ = null, sinσ = null, cosσ = null; // σ = angular distance P₁ P₂ on the sphere
let cos2σₘ = null;                      // σₘ = angular distance on the sphere from the equator to the midpoint of the line
let cosSqα = null;                      // α = azimuth of the geodesic at the equator

let λʹ = null;
do {
sinλ = Math.sin(λ);
cosλ = Math.cos(λ);
const sinSqσ = (cosU2*sinλ) * (cosU2*sinλ) + (cosU1*sinU2-sinU1*cosU2*cosλ)**2);
sinσ = Math.sqrt(sinSqσ);
cosσ = sinU1*sinU2 + cosU1*cosU2*cosλ;
σ = Math.atan2(sinσ, cosσ);
const sinα = cosU1 * cosU2 * sinλ / sinσ;
cosSqα = 1 - sinα*sinα;
cos2σₘ = cosσ - 2*sinU1*sinU2/cosSqα;
const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
λʹ = λ;
λ = L + (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
} while (Math.abs(λ-λʹ) > 1e-12);

const uSq = cosSqα * (a*a - b*b) / (b*b);
const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));

const s = b*A*(σ-Δσ); // s = length of the geodesic

const α1 = Math.atan2(cosU2*sinλ,  cosU1*sinU2-sinU1*cosU2*cosλ); // initial bearing
const α2 = Math.atan2(cosU1*sinλ, -sinU1*cosU2+cosU1*sinU2*cosλ); // final bearing
```

See below for the complete code.

### Destination given distance & bearing from start point (direct solution)

Again, this presents Vincenty’s formulæ arranged to be close to how they are used in the script.

 a, b = major & minor semi-axes of the ellipsoid f = flattening (a−b)/a φ1, φ2 = geodetic latitude L = difference in longitude s = length of the geodesic along the surface of the ellipsoid (in the same units as a & b) α1, α2 = azimuths of the geodesic (initial/final bearing) tan U1 = (1−f) · tan φ1 (U is ‘reduced latitude’) cos U1 = 1 / √1 + tan² U1, sin U1 = tan U1 · cos U1 (trig identities; §6) σ1 = atan(tan U1 / cos α1) (1) sin α = cos U1 · sin α1 (2) cos² α = 1 − sin² α (trig identity; §6) 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) σ = s / (b·A) (first approximation) iterate until change in σ is negligible (e.g. 10-12 ≈ 0.006mm) { cos 2σm = cos(2σ1 + σ) (5) Δσ = B · sin σ · {cos 2σm + B/4 · [cos σ · (−1 + 2 · cos² 2σm)      − B/6 · cos 2σm · (−3 + 4 · sin² σ) · (−3 + 4 · cos² 2σm)]} (6) σʹ = s / b·A + Δσ (7) } φ2 = atan(sin U1 · cos σ + cos U1 · sin σ · cos α1      / (1−f) · √sin² α + (sin U1 · sin σ − cos U1 · cos σ · cos α1)² ) (8) λ = atan(sin σ · sin α1 / cos U1 · cos σ − sin U1 · sin σ · cos α1) (9) C = f/16 · cos² α · [4 + f · (4 − 3 · cos² α)] (10) L = λ − (1−C) · f · sin α · {σ + C · sin σ · [cos 2σm + C · cos σ · (−1 + 2 · cos² 2σm)]} (11) λ2 = λ1 + L α2 = atan( sin α / −(sin U1 · sin σ − cos U1 · cos σ · cos α1) ) (12)

Where:

• φ2, λ2 is destination point
• α2 is final bearing (in direction p1→p2)

#### JavaScript

(slightly simplified, ignoring convergence etc)
```const sinα1 = Math.sin(α1);
const cosα1 = Math.cos(α1);

const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
const σ1 = Math.atan2(tanU1, cosα1); // σ1 = angular distance on the sphere from the equator to P1
const sinα = cosU1 * sinα1;          // α = azimuth of the geodesic at the equator
const cosSqα = 1 - sinα*sinα;
const uSq = cosSqα * (a*a - b*b) / (b*b);
const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

let σ = s / (b*A), sinσ = null, cosσ = null; // σ = angular distance P₁ P₂ on the sphere
let cos2σₘ = null; // σₘ = angular distance on the sphere from the equator to the midpoint of the line

let σʹ = null;
do {
cos2σₘ = Math.cos(2*σ1 + σ);
sinσ = Math.sin(σ);
cosσ = Math.cos(σ);
const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));
σʹ = σ;
σ = s / (b*A) + Δσ;
} while (Math.abs(σ-σʹ) > 1e-12);

const x = sinU1*sinσ - cosU1*cosσ*cosα1;
const φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + x*x));
const λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1);
const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
const λ2 = λ1 + L;

const α2 = Math.atan2(sinα, -x); // final bearing
```

See below for the complete code.

#### Datums

While is it wonderful to be able to calculate distances to millimetre accuracy, there are certain caveats which it is worth being aware of if you are looking for accuracy better than around ±1 metre.

Millimetre precision is stretching the limits of geodetic survey techniques: very few points will have been surveyed to such precision (it’s easy to provide lati­tude/longi­tude ref­erences to 6 decimals of degrees; less easy to make that precision meaningful rather than spurious).

Further, the familiar WGS-84 (World Geodetic System 1984) has no ‘physical realisa­tion’ – it is not tied to geodetic groundsta­tions, just to satellites – and is defined to be accurate to no better than around ±1 metre, irrespective of how good your GPS unit is [actually, the ‘G-series’ are a little more accurate, but there are no official transforma­tion param­eters]. A metre is good enough for most purposes, but a long way from Vincenty’s millimetre accuracy.

Then, when working to an accuracy of better than ±1 metre, plate tectonic movements become significant. Depending on the datum it’s defined in, any given point will most likely have shifted many millimetres from where it was last year.

Simplifying somewhat (well, actually, a lot!), in order to address requirements for greater precision, the ITRS was developed to give greater accuracy than WGS84, with ‘epoch’-dependant ITRF ‘realisa­tions’, where the lati­tude/longi­tude coordinate of a posi­tion will vary over time.

Various ‘static’ reference frames are also defined for various continents – NAD83 for North America, ETRS89 for Europe, GDA94 for Australia, etc – within which lati­tude/longi­tude coordinates remain essentially fixed (at least to centimetre or so accuracy, major earthquake events excepted). These reference frames have ‘epoch’-dependant mappings to ITRF datums. (Due to plate tectonic continental drift, ETRS89 shifts against ITRF by about 25mm/year; GDA94 by around 80mm/year).

And for any latitude/longitude coordinates dating from earlier datums (such as NAD27, ED50, NTF, AGD66, and any number of others), even metre precision is exceeding their limits.

Also, it must be remembered that these distances are calculated on the surface of the reference ellipsoid, which will differ (to varying degrees) from the earth geoid, and even more from the earth’s physical surface. The geoid is equivalent to the local mean-sea-level (or its equivalent in land-locked locales), and is determined by the local force of gravity; it varies in height from the ellipsoid by around ±100 metres. The earth’s surface varies much more; if you happen to be located in Colorado, 2 kilometres above msl, distances will be 0.03% greater than on the ellipsoid surface; +1 metre over just 30 kilometres.

So for millimetre accuracy to be relevant, you need to know a fair bit about the reference frames, datums, and epochs of your coordinates, and other param­eters of your measure­ments. Caveat emptor!

The ubiquitous WGS-84 is a geocentric datum, based on an ellipsoid with:

 Semi-major axis a = 6 378 137.0 metres Semi-minor axis b ≈ 6 356 752.314 245 metres Inverse flattening 1/f = 298.257 223 563

Latitude/longitude points can be converted between different datums by converting them to geocentric cartesian coordinates, applying a 7-param­eter Helmert transforma­tion, and converting back to a geodetic coordinate in the new datum.

Should even millimetre accuracy not be sufficient for you, Charles Karney has improved on Vincenty’s method with a method which has errors in nanometers (and always converges on antipodal points) – Geodesics on an ellipsoid of revolu­tion, 2011, Algorithms for geodesics, 2012; though Karney’s method is a lot less straight­forward (someday I may get around to an implementa­tion).

### Live examples

... to verify the code is correct!

Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:

• deg-min-sec suffixed with N/S/E/W (e.g. 40°44′55″N, 73 59 11W), or
• signed decimal degrees without compass direction, where negative indicates west/south (e.g. 40.7486, -73.9864):
Inverse solution
 Point 1: , Point 2: ,
 Distance: m Initial bearing: Final bearing:
Direct solution
 Start point: , Bearing: Distance: m
 Destination point: Final bearing:

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)

• 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 = rad⋅180/π), 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 code involving atan2 in VBA (using Excel.WorksheetFunction.Atan2), you will need to reverse the arguments, as Excel has them the opposite way around from JavaScript – conven­tional order is atan2(y, x), but Excel uses atan2(x, y). To use this code in VBA, you will also have to rename the variables A, B / a, b, as VBA is case-insensitive (and replace greek letters with Latin equivalents?).
• 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 JavaScript source code, also available on GitHub. Full documentation is available, as well as a test suite.

Note I use Greek letters in variables representing maths symbols conven­tionally presented as Greek letters: I value the great benefit in legibility over the minor inconvenience in typing (if you encounter any problems, ensure your `<head>` includes `<meta charset="utf-8">`), and/or use UTF-8 encoding when saving files.

With its untyped C-style syntax, JavaScript reads remarkably close to pseudo-code: exposing the algorithms with a minimum of syntactic distrac­tions. These func­tions should be simple to translate into other languages if required, though can also be used as-is in browsers and Node.js.

For convenience & clarity, I have extended the base JavaScript `Number` object with `toRadians()` and `toDegrees()` methods: I don’t see great likelihood of conflicts, as these are ubiquitous operations.

February 2019: I have refactored the library to use ES modules, as well as extending it in scope; see the GitHub README and CHANGELOG for details.

Other languages: I cannot support translations in languages I’m not familiar with, but if you have ported the code to another language, I am happy to provide links here.

• Jean Brouwers has made a Python version.
• Tomasz Jastrzębski has implemented it for Excel. I offer these scripts for free use and adapta­tion to balance my debt to the open-source info-verse. You are welcome to re-use these scripts [under an MIT licence, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

If you need any advice or development work done, I am available for consultancy.

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