Fix d3.geo.interpolate for coincident points.

This also improves accuracy for small distances, using the haversine
distance formula instead of the law of cosines.

Fixes d3#1080.

Makefile

@@ -198,6 +198,7 @@ d3.geo.js: \
 	src/geo/equirectangular.js \
 	src/geo/gnomonic.js \
 	src/geo/graticule.js \
+	src/geo/haversin.js \
 	src/geo/interpolate.js \
 	src/geo/greatArc.js \
 	src/geo/mercator.js \

d3.js

@@ -6130,15 +6130,20 @@ d3 = function() {
       });
     };
   }
+  function d3_geo_haversin(x) {
+    return (x = Math.sin(x / 2)) * x;
+  }
   d3.geo.interpolate = function(source, target) {
     return d3_geo_interpolate(source[0] * d3_radians, source[1] * d3_radians, target[0] * d3_radians, target[1] * d3_radians);
   };
   function d3_geo_interpolate(x0, y0, x1, y1) {
-    var cy0 = Math.cos(y0), sy0 = Math.sin(y0), cy1 = Math.cos(y1), sy1 = Math.sin(y1), kx0 = cy0 * Math.cos(x0), ky0 = cy0 * Math.sin(x0), kx1 = cy1 * Math.cos(x1), ky1 = cy1 * Math.sin(x1), d = Math.acos(Math.max(-1, Math.min(1, sy0 * sy1 + cy0 * cy1 * Math.cos(x1 - x0)))), k = 1 / Math.sin(d);
-    function interpolate(t) {
+    var cy0 = Math.cos(y0), sy0 = Math.sin(y0), cy1 = Math.cos(y1), sy1 = Math.sin(y1), kx0 = cy0 * Math.cos(x0), ky0 = cy0 * Math.sin(x0), kx1 = cy1 * Math.cos(x1), ky1 = cy1 * Math.sin(x1), d = 2 * Math.asin(Math.sqrt(d3_geo_haversin(y1 - y0) + cy0 * cy1 * d3_geo_haversin(x1 - x0))), k = 1 / Math.sin(d);
+    var interpolate = d ? function(t) {
       var B = Math.sin(t *= d) * k, A = Math.sin(d - t) * k, x = A * kx0 + B * kx1, y = A * ky0 + B * ky1, z = A * sy0 + B * sy1;
-      return [ Math.atan2(y, x) / d3_radians, Math.atan2(z, Math.sqrt(x * x + y * y)) / d3_radians ];
-    }
+      return [ Math.atan2(y, x) * d3_degrees, Math.atan2(z, Math.sqrt(x * x + y * y)) * d3_degrees ];
+    } : function() {
+      return [ x0 * d3_degrees, y0 * d3_degrees ];
+    };
     interpolate.distance = d;
     return interpolate;
   }

src/geo/haversin.js

@@ -0,0 +1,3 @@
+function d3_geo_haversin(x) {
+  return (x = Math.sin(x / 2)) * x;
+}

src/geo/interpolate.js

@@ -14,20 +14,20 @@ function d3_geo_interpolate(x0, y0, x1, y1) {
       ky0 = cy0 * Math.sin(x0),
       kx1 = cy1 * Math.cos(x1),
       ky1 = cy1 * Math.sin(x1),
-      d = Math.acos(Math.max(-1, Math.min(1, sy0 * sy1 + cy0 * cy1 * Math.cos(x1 - x0)))),
+      d = 2 * Math.asin(Math.sqrt(d3_geo_haversin(y1 - y0) + cy0 * cy1 * d3_geo_haversin(x1 - x0))),
       k = 1 / Math.sin(d);
 
-  function interpolate(t) {
+  var interpolate = d ? function(t) {
     var B = Math.sin(t *= d) * k,
         A = Math.sin(d - t) * k,
         x = A * kx0 + B * kx1,
         y = A * ky0 + B * ky1,
         z = A * sy0 + B * sy1;
     return [
-      Math.atan2(y, x) / d3_radians,
-      Math.atan2(z, Math.sqrt(x * x + y * y)) / d3_radians
+      Math.atan2(y, x) * d3_degrees,
+      Math.atan2(z, Math.sqrt(x * x + y * y)) * d3_degrees
     ];
-  }
+  } : function() { return [x0 * d3_degrees, y0 * d3_degrees]; };
 
   interpolate.distance = d;
 

test/geo/interpolate-test.js

@@ -0,0 +1,25 @@
+require("../env");
+
+var vows = require("vows"),
+    assert = require("assert");
+
+var suite = vows.describe("d3.geo.interpolate");
+
+suite.addBatch({
+  "interpolate": {
+    topic: function() {
+      return d3.geo.interpolate;
+    },
+    "zero distance": function(arc) {
+      assert.deepEqual(d3.geo.interpolate([140.63289, -29.95101], [140.63289, -29.95101])(.5), [140.63289, -29.95101]);
+    },
+    "equator": function(arc) {
+      assert.inDelta(d3.geo.interpolate([10, 0], [20, 0])(.5), [15, 0], 1e-6);
+    },
+    "meridian": function(arc) {
+      assert.inDelta(d3.geo.interpolate([10, -20], [10, 40])(.5), [10, 10], 1e-6);
+    }
+  }
+});
+
+suite.export(module);