--- a/flys-artifacts/ChangeLog	Wed Apr 20 13:13:19 2011 +0000
+++ b/flys-artifacts/ChangeLog	Wed Apr 20 14:26:51 2011 +0000
@@ -1,3 +1,17 @@
+2011-04-20	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java:
+	  Implementation of "Abflusskurve/Abflusstafel" calculation.
+
+	  Added method interpolateWQ() which takes an km and results in a
+	  tuple of two double arrays containing the w/q values interpolated
+	  between the surrounding w/q values of the table.
+	  w values are interpolated linear, q values with a cubic spline.
+
+	  Drawing w over q gives you the discharge table at the given km.
+
+	  !!! This code needs testing !!!
+
 2011-04-20	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* pom.xml: Added dependency to Apache Commons Math 2.2 (Apache License 2.0)
@@ -15,7 +29,7 @@
 2011-04-19	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java:
-	  add a method interpolateW() method which takes an array of
+	  add a method interpolateW() which takes an array of
 	  q values and returns an equal sized array of w values.
 	  This is essentially the "Wasserstand/Wasserspiegellagen" calculation
 	  of desktop FLYS.

--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Wed Apr 20 13:13:19 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Wed Apr 20 14:26:51 2011 +0000
@@ -15,6 +15,14 @@
 import java.util.Collections;
 import java.util.Iterator;
 
+import org.apache.log4j.Logger;
+
+import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
+
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+import org.apache.commons.math.ArgumentOutsideDomainException;
+
 import org.hibernate.Session;
 import org.hibernate.Query;
 import org.hibernate.SQLQuery;
@@ -24,6 +32,7 @@
 public class WstValueTable
 implements   Serializable
 {
+    private static Logger log = Logger.getLogger(WstValueTable.class);
 
     // TODO: put this into a property file
     public static final String SQL_POS_WQ = 
@@ -31,6 +40,8 @@
         "    FROM wst_value_table"          +
         "    WHERE wst_id = :wst_id";
 
+    public static final double DEFAULT_STEP_WIDTH = 0.01;
+
     public static class Column
     implements          Serializable
     {
@@ -188,6 +199,134 @@
             if (d > 0.0) return +1;
             return 0;
         }
+
+        public double [][] cloneWQs() {
+            return new double [][] {
+                (double [])ws.clone(),
+                (double [])ws.clone() };
+        }
+
+        public double [][] interpolateWQ(Row other, double km, double stepWidth) {
+
+            int W1 =       ascendingWs();
+            int W2 = other.ascendingWs();
+
+            int W = Math.min(W1, W2);
+
+            if (W < 1) {
+                return new double[2][0];
+            }
+
+            double factor = factor(km, this.km, other.km);
+
+            double minW = weight(factor, ws[0],   other.ws[0]);
+            double maxW = weight(factor, ws[W-1], other.ws[W-1]);
+
+            if (minW > maxW) {
+                double t = minW;
+                minW = maxW;
+                maxW = t;
+            }
+            double [] x = new double[W];
+            double [] y = new double[W];
+
+            for (int i = 0; i < W; ++i) {
+                x[i] = weight(factor, ws[i], other.ws[i]);
+                y[i] = weight(factor, qs[i], other.qs[i]);
+            }
+
+            SplineInterpolator interpolator = new SplineInterpolator();
+            PolynomialSplineFunction spline = interpolator.interpolate(x, y);
+
+            double [] outWs =
+                new double[(int)Math.ceil((maxW - minW)/stepWidth)];
+            double [] outQs =
+                new double[outWs.length];
+
+            try {
+                double w = minW;
+                for (int i = 0; i < outWs.length; ++i, w += stepWidth) {
+                    outQs[i] = spline.value(outWs[i] = w);
+                }
+            }
+            catch (ArgumentOutsideDomainException aode) {
+                log.error("Spline interpolation failed.", aode);
+            }
+
+            return new double [][] { outWs, outQs };
+
+        }
+
+        public double [][] interpolateWQ(double stepWidth) {
+            int W = ascendingWs(); // ignore back jumps
+
+            if (W < 1) {
+                return new double[2][0];
+            }
+
+            double [] x = new double[W];
+            double [] y = new double[W];
+
+            for (int i = 0; i < W; ++i) {
+                x[i] = ws[i];
+                y[i] = qs[i];
+            }
+
+            SplineInterpolator interpolator = new SplineInterpolator();
+
+            PolynomialSplineFunction spline = interpolator.interpolate(x, y);
+
+            double minW = ws[0];
+            double maxW = ws[W-1];
+
+            double [] outWs =
+                new double[(int)Math.ceil((maxW - minW)/stepWidth)];
+            double [] outQs =
+                new double[outWs.length];
+
+            try {
+                double w = minW;
+                for (int i = 0; i < outWs.length; ++i, w += stepWidth) {
+                    outQs[i] = spline.value(outWs[i] = w);
+                }
+            }
+            catch (ArgumentOutsideDomainException aode) {
+                log.error("Spline interpolation failed.", aode);
+            }
+
+            return new double [][] { outWs, outQs };
+        }
+
+        public int ascendingWs() {
+            if (ws.length < 2) {
+                return ws.length;
+            }
+
+            int idx = 1;
+
+            for (; idx < ws.length; ++idx) {
+                if (Double.isNaN(ws[idx]) || ws[idx] < ws[idx-1]) {
+                    return idx;
+                }
+            }
+
+            return idx;
+        }
+
+        public double [][] weightWQs(Row other, double km) {
+            int W = Math.min(ws.length, other.ws.length);
+            double factor = factor(km, this.km, other.km);
+
+            double [] outWs = new double[W];
+            double [] outQs = new double[W];
+
+            for (int i = 0; i < W; ++i) {
+                outWs[i] = weight(factor, ws[i], other.ws[i]);
+                outQs[i] = weight(factor, qs[i], other.qs[i]);
+            }
+
+            return new double [][] { outWs, outQs };
+        }
     }
     // class Row
 
@@ -234,6 +373,40 @@
         return ws;
     }
 
+
+    public double [][] interpolateWQ(double km) {
+        return interpolateWQ(km, DEFAULT_STEP_WIDTH, true);
+    }
+
+    public double [][] interpolateWQ(
+        double  km,
+        double  stepWidth, 
+        boolean checkAscending
+    ) {
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        if (rowIndex >= 0) { // direct row match
+            Row row = rows.get(rowIndex);
+            return checkAscending
+                ? row.interpolateWQ(stepWidth)
+                : row.cloneWQs();
+        }
+
+        rowIndex = -rowIndex -1;
+
+        if (rowIndex < 1 || rowIndex >= rows.size()) {
+            // do not extrapolate
+            return new double[2][0];
+        }
+
+        Row r1 = rows.get(rowIndex-1);
+        Row r2 = rows.get(rowIndex);
+
+        return checkAscending
+            ? r1.interpolateWQ(r2, km, stepWidth)
+            : r1.weightWQs(r2, km);
+    }
+
     public static WstValueTable getTable(River river) {
         return getTable(river, 0);
     }