From bacb0a917cf3dc9817a5c9554a59d73be2413fc7 Mon Sep 17 00:00:00 2001
From: crschmidt <crschmidt@crschmidt.net>
Date: Thu, 8 Jun 2006 03:37:20 +0000
Subject: [PATCH] Add vincenty great circle formula.

git-svn-id: http://svn.openlayers.org/trunk/openlayers@554 dc9f47b5-9b13-0410-9fdd-eb0c1a62fdaf
---
 lib/OpenLayers/Util.js | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/lib/OpenLayers/Util.js b/lib/OpenLayers/Util.js
index f7f1211e0a..1f5420b657 100644
--- a/lib/OpenLayers/Util.js
+++ b/lib/OpenLayers/Util.js
@@ -949,3 +949,39 @@ OpenLayers.Util.mouseLeft = function (evt, div) {
     // if the target we stop at isn't the div, then we've left the div.
     return (target != div);
 };
+
+OpenLayers.Util.rad = function(x) {return x*Math.PI/180;};
+OpenLayers.Util.distVincenty=function(p1, p2) {
+    var a = 6378137, b = 6356752.3142,  f = 1/298.257223563;
+    var L = OpenLayers.Util.rad(p2.lon - p1.lon);
+    var U1 = Math.atan((1-f) * Math.tan(OpenLayers.Util.rad(p1.lat)));
+    var U2 = Math.atan((1-f) * Math.tan(OpenLayers.Util.rad(p2.lat)));
+    var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
+    var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
+    var lambda = L, lambdaP = 2*Math.PI;
+    var iterLimit = 20;
+    while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0) {
+        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 alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma);
+        var cosSqAlpha = Math.cos(alpha) * Math.cos(alpha);
+        var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
+        var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
+        lambdaP = lambda;
+        lambda = L + (1-C) * f * Math.sin(alpha) *
+        (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
+    }
+    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);
+    var d = s.toFixed(3)/1000; // round to 1mm precision
+    return d;
+};