Vincenty ‘Direct’ formula | Movable Type Scripts

In addition to the ‘inverse’ formula for calculating distance between two points on an accurate ellipsoidal model of the earth, Thaddeus Vincenty devised a formula for deriving the destination point given a start point, an initial bearing, and a distance travelled. This is an implementation of that formula in JavaScript. For more details, see the ‘inverse’ page with the script for determining the distance between two points.

Enter the co-ordinates into the text boxes to try it out (using deg-min-sec suffixed with N/S/E/W, or signed decimal degrees):

Start Lat:		Start Long:
Bearing:		Distance:

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

a, b = major & minor semiaxes of the ellipsoid
f = flattening (a−b)/a
φ₁, φ₂ = geodetic latitude
s = length of the geodesic
α₁, α₂ = azimuths of the geodesic (initial/final bearing)

tanU₁ = (1−f).tanφ₁ (U is ‘reduced latitude’)
cosU₁ = 1/√(1+tan²U₁), sinU₁ = tanU₁.cosU₁ (trig identities; §6)
σ₁ = atan2(tanU₁, cosα₁)	(1)
sinα = cosU₁.sinα₁	(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 (1st approximation), σ′ = 2π
while abs(σ−σ′) > 10^-12 { (i.e. 0.06mm)
cos2σ_m = cos(2.σ₁ + σ)	(5)
Δσ = 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 + Δσ	(7)
}
φ₂ = atan2(sinU₁.cosσ + cosU₁.sinσ.cosα₁, (1−f).√[sin²α + (sinU₁.sinσ − cosU₁.cosσ.cosα₁)²])	(8)
λ = atan2(sinσ.sinα₁, cosU₁.cosσ − sinU₁.sinσ.cosα₁)	(9)
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)]} (difference in longitude)	(11)
α₂ = atan2(sinα, −sinU₁.sinσ + cosU₁.cosσ.cosα₁) (reverse azimuth)	(12)
p₂ = (φ₂, λ₁+L)

See below for the source code of the JavaScript implementation. The JavaScript should be quite straightforward to translate into other languages if required (perhaps with alternate approach to returning latitude/longitude object).

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.


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Vincenty Direct Solution of Geodesics on the Ellipsoid (c) Chris Veness 2005-2012              */
/*                                                                                                */
/* from: Vincenty direct 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 destination point given start point lat/long, bearing & distance, 
 * using Vincenty inverse formula for ellipsoids
 *
 * @param   {Number} lat1, lon1: first point in decimal degrees
 * @param   {Number} brng: initial bearing in decimal degrees
 * @param   {Number} dist: distance along bearing in metres
 * @returns (LatLon} destination point
 */
function destVincenty(lat1, lon1, brng, dist) {
  var a = 6378137, b = 6356752.3142,  f = 1/298.257223563;  // WGS-84 ellipsiod
  var s = dist;
  var alpha1 = brng.toRad();
  var sinAlpha1 = Math.sin(alpha1);
  var cosAlpha1 = Math.cos(alpha1);
  
  var tanU1 = (1-f) * Math.tan(lat1.toRad());
  var cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1*cosU1;
  var sigma1 = Math.atan2(tanU1, cosAlpha1);
  var sinAlpha = cosU1 * sinAlpha1;
  var cosSqAlpha = 1 - sinAlpha*sinAlpha;
  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 sigma = s / (b*A), sigmaP = 2*Math.PI;
  while (Math.abs(sigma-sigmaP) > 1e-12) {
    var cos2SigmaM = Math.cos(2*sigma1 + sigma);
    var sinSigma = Math.sin(sigma);
    var cosSigma = Math.cos(sigma);
    var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
      B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
    sigmaP = sigma;
    sigma = s / (b*A) + deltaSigma;
  }

  var tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1;
  var lat2 = Math.atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1, 
      (1-f)*Math.sqrt(sinAlpha*sinAlpha + tmp*tmp));
  var lambda = Math.atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1);
  var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
  var L = lambda - (1-C) * f * sinAlpha *
      (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
  var lon2 = (lon1.toRad()+L+3*Math.PI)%(2*Math.PI) - Math.PI;  // normalise to -180...+180

  var revAz = Math.atan2(sinAlpha, -tmp);  // final bearing, if required

  return { lat: lat2.toDeg(), lon: lon2.toDeg(), finalBearing: revAz.toDeg() };
}

// ---- extend Number object with methods for converting degrees/radians

Number.prototype.toRad = function() {
  return this * Math.PI / 180;
}

Number.prototype.toDeg = function() {
  return this * 180 / Math.PI;
}

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