Makes splineMonotone() actually monotone; re: #91

Author: danburzo

src/interpolate/splineBasis.js

@@ -1,7 +1,6 @@
 import identity from '../util/identity';
 
 /*
-
 	Basis spline
 	------------
 

src/interpolate/splineMonotone.js

@@ -1,19 +1,46 @@
 import identity from '../util/identity';
 
-const monotone = (v0, v1, v2, v3, h, t) => {
+/* 
+	Monotone spline
+	---------------
+
+	Based on:
+
+		Steffen, M.
+		"A simple method for monotonic interpolation in one dimension."
+		in Astronomy and Astrophysics, Vol. 239, p. 443-450 (Nov. 1990),
+      	Provided by the SAO/NASA Astrophysics Data System.
+
+		https://ui.adsabs.harvard.edu/abs/1990A&A...239..443S
+
+	(Reference thanks to `d3/d3-shape`)
+*/
+
+const sgn = Math.sign,
+	min = Math.min,
+	abs = Math.abs;
+
+const monotone = (y_im1, y_i, y_ip1, y_ip2, h, t) => {
 	let h2 = h * h;
 	let t2 = t * t;
 	let t3 = t2 * t;
 
-	let s20 = (v2 - v0) / (2 * h);
-	let s31 = (v3 - v1) / (2 * h);
-	let s21 = (v2 - v1) / h;
+	let s_im1 = (y_i - y_im1) / h;
+	let s_i = (y_ip1 - y_i) / h;
+	let s_ip1 = (y_ip2 - y_ip1) / h;
+
+	let yp_i =
+		(sgn(s_im1) + sgn(s_i)) *
+		min(abs(s_im1), abs(s_i), 0.5 * abs((y_ip1 - y_im1) / (2 * h)));
+	let yp_ip1 =
+		(sgn(s_i) + sgn(s_ip1)) *
+		min(abs(s_i), abs(s_ip1), 0.5 * abs((y_ip2 - y_i) / (2 * h)));
 
 	return (
-		((s20 + s31 - 2 * s21) / h2) * t3 +
-		((3 * s21 - 2 * s20 - s31) / h) * t2 +
-		s20 * t +
-		v1
+		((yp_i + yp_ip1 - 2 * si) / h2) * t3 +
+		((3 * si - 2 * yp_i - yp_ip1) / h) * t2 +
+		yp_i * t +
+		y_i
 	);
 };