Vincenty ‘Direct’ formula
 

 

In addition to the 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:

The JavaScript should be quite straightforward to translate into other languages, if required.

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  
s = length of the geodesic  
α1, α2 = azimuths of the geodesic (initial/final bearing)  
   
tanU1 = (1−f).tanφ1 (U is ‘reduced latitude’)  
cosU1 = 1/√(1+tan²U1), sinU1 = tanU1.cosU1 (trig identities; §6)  
σ1 = atan2(tanU1, cosα1) (1)
sinα = cosU1.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 (1st approximation), σ′ = 2π  
while abs(σ−σ′) > 10-12 { (i.e. 0.06mm)  
     cos2σm = cos(2.σ1 + σ) (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)
}    
φ2 = atan2(sinU1.cosσ + cosU1.sinσ.cosα1, (1−f).√[sin²α + (sinU1.sinσ − cosU1.cosσ.cosα1)²]) (8)
λ = atan2(sinσ.sinα1, cosU1.cosσ − sinU1.sinσ.cosα1) (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)
α2 = atan(sinα, −sinU1.sinσ + cosU1.cosσ.cosα1) (reverse azimuth) (12)
p2 = (φ2, λ1+L)  
   
  JavaScript function is shown below. View Source to see complete implementation. You are welcome to re-use these scripts [without any warranty express or implied] provided you retain my copyright notice and when possible a link to my website. If you have any queries or find any problems, please contact me.
© 2002-2006 Chris Veness
    

/*
 * Calculate destination point given start point lat/long (numeric degrees), 
 * bearing (numeric degrees) & distance (in m).
 *
 * 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
 */
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), 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), 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 revAz = Math.atan2(sinAlpha, -tmp);  // final bearing

  return new LatLon(lat2.toDeg(), lon1+L.toDeg());
}