From 49a53f19358a5f6233a4925f6a782407a1d3d7fe Mon Sep 17 00:00:00 2001
From: Xavier Raynaud <xavier.raynaud@sintef.no>
Date: Thu, 15 Mar 2012 16:15:32 +0100
Subject: [PATCH] Compute the derivative directly (not using epsilon
 parameter).

---
 opm/core/utility/linearInterpolation.hpp | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/opm/core/utility/linearInterpolation.hpp b/opm/core/utility/linearInterpolation.hpp
index 7b40a5c8..0ed8b57b 100644
--- a/opm/core/utility/linearInterpolation.hpp
+++ b/opm/core/utility/linearInterpolation.hpp
@@ -74,12 +74,15 @@ namespace Opm
 				    const std::vector<T>& yv,
 				    double x)
     {
-	double epsilon = 1e-4; // @@ Ad hoc, should choose based on xv.
-	double x_low = std::max(xv[0], x - epsilon);
-	double x_high = std::min(xv.back(), x + epsilon);
-	T low = linearInterpolation(xv, yv, x_low);
-	T high = linearInterpolation(xv, yv, x_high);
-	return (high - low)/(x_high - x_low);
+	std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
+	int lb_ix = lb - xv.begin();
+	if (lb_ix == 0) {
+	    return 0.;
+	} else if (lb_ix == int(xv.size())) {
+	    return 0.;
+	} else {
+	    return (yv[lb_ix] - yv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
+	}
     }