Changeset 554

Timestamp:

06/07/06 23:37:20 (2 years ago)

Author:

crschmidt

Message:

Add vincenty great circle formula.

Files:

• trunk/openlayers/lib/OpenLayers/Util.js (modified) (1 diff)

trunk/openlayers/lib/OpenLayers/Util.js

@@ -950 +950 @@
-    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;
+};