From 60e869e295c67e714fd26bca4e39757d85185de7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=89ric=20Lemoine?= <eric.lemoine@camptocamp.com>
Date: Thu, 24 Jun 2010 09:05:36 +0000
Subject: [PATCH] Add Vincenty direct formula, and make the Vincenty functions
 part of the API, p=aabt, r=me (closes #2692)

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

diff --git a/lib/OpenLayers/Util.js b/lib/OpenLayers/Util.js
index 8961dd5c65..030ad4ae0f 100644
--- a/lib/OpenLayers/Util.js
+++ b/lib/OpenLayers/Util.js
@@ -892,22 +892,45 @@ OpenLayers.Util.toFloat = function (number, precision) {
 OpenLayers.Util.rad = function(x) {return x*Math.PI/180;};
 
 /**
- * Function: distVincenty
+ * Function: deg
+ *
+ * Parameters:
+ * x - {Float}
+ *
+ * Returns:
+ * {Float}
+ */
+OpenLayers.Util.deg = function(x) {return x*180/Math.PI;};
+
+/**
+ * Property: VincentyConstants
+ * {Object} Constants for Vincenty functions.
+ */
+OpenLayers.Util.VincentyConstants = {
+    a: 6378137,
+    b: 6356752.3142,
+    f: 1/298.257223563
+};
+
+/**
+ * APIFunction: distVincenty
  * Given two objects representing points with geographic coordinates, this
  *     calculates the distance between those points on the surface of an
  *     ellipsoid.
- * 
+ *
  * Parameters:
  * p1 - {<OpenLayers.LonLat>} (or any object with both .lat, .lon properties)
  * p2 - {<OpenLayers.LonLat>} (or any object with both .lat, .lon properties)
- * 
+ *
  * Returns:
  * {Float} The distance (in km) between the two input points as measured on an
  *     ellipsoid.  Note that the input point objects must be in geographic
  *     coordinates (decimal degrees) and the return distance is in kilometers.
  */
-OpenLayers.Util.distVincenty=function(p1, p2) {
-    var a = 6378137, b = 6356752.3142,  f = 1/298.257223563;
+OpenLayers.Util.distVincenty = function(p1, p2) {
+    var ct = OpenLayers.Util.VincentyConstants;
+    var a = ct.a, b = ct.b, f = ct.f;
+
     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)));
@@ -945,6 +968,68 @@ OpenLayers.Util.distVincenty=function(p1, p2) {
     return d;
 };
 
+/**
+ * APIFunction: destinationVincenty
+ * Calculate destination point given start point lat/long (numeric degrees),
+ * bearing (numeric degrees) & distance (in m).
+ * Adapted from Chris Veness work, see
+ * http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html
+ *
+ * Parameters:
+ * lonlat  - {<OpenLayers.LonLat>} (or any object with both .lat, .lon
+ *     properties) The start point.
+ * brng     - {Float} The bearing (degrees).
+ * distance - {Float} The ground distance (meters).
+ *
+ * Returns:
+ * {<OpenLayers.LonLat>} The destination point.
+ */
+OpenLayers.Util.destinationVincenty = function(lonlat, brng, dist) {
+    var u = OpenLayers.Util;
+    var ct = u.VincentyConstants;
+    var a = ct.a, b = ct.b, f = ct.f;
+
+    var lon1 = lonlat.lon;
+    var lat1 = lonlat.lat;
+
+    var s = dist;
+    var alpha1 = u.rad(brng);
+    var sinAlpha1 = Math.sin(alpha1);
+    var cosAlpha1 = Math.cos(alpha1);
+
+    var tanU1 = (1-f) * Math.tan(u.rad(lat1));
+    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 revAz = Math.atan2(sinAlpha, -tmp);  // final bearing
+
+    return new OpenLayers.LonLat(lon1+u.deg(L), u.deg(lat2));
+};
+
 /**
  * Function: getParameters
  * Parse the parameters from a URL or from the current page itself into a