--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java	Tue Feb 14 09:02:51 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java	Tue Feb 14 16:48:13 2012 +0000
@@ -4,7 +4,6 @@
 import java.util.Map;
 import java.util.HashMap;
 import java.util.Arrays;
-import java.util.Comparator;
 
 import java.io.Serializable;
 
@@ -18,6 +17,8 @@
 import de.intevation.flys.model.DischargeTable;
 import de.intevation.flys.model.DischargeTableValue;
 
+import gnu.trove.TDoubleArrayList;
+
 public class DischargeTables
 implements   Serializable
 {
@@ -124,42 +125,9 @@
             // TODO: Filter by time interval
             DischargeTable table = tables.get(0);
 
-            boolean qSorted = true;
-
-            final double[][] vs =
-                loadDischargeTableValues(table, scale, qSorted);
-
-            if (!qSorted) {
-                log.debug("need to sort by q values");
-                // TODO: Do this db level.
-                // XXX: This is so ugly :-(
-                Integer [] indices = new Integer[vs[0].length];
-                for (int i = 0; i < indices.length; ++i) {
-                    indices[i] = i;
-                }
+            double [][] vs = loadDischargeTableValues(table, scale);
 
-                Arrays.sort(indices, new Comparator<Integer>() {
-                    public int compare(Integer a, Integer b) {
-                        double va = vs[1][a];
-                        double vb = vs[1][b];
-                        double d = va - vb;
-                        if (d < 0.0) return -1;
-                        if (d > 0.0) return +1;
-                        return 0;
-                    }
-                });
-
-                double [][] vs2 = new double[2][vs[0].length];
-                for (int i = 0; i < indices.length; ++i) {
-                    vs2[0][i] = vs[0][indices[i]];
-                    vs2[1][i] = vs[1][indices[i]];
-                }
-                values.put(gaugeName, vs2);
-            }
-            else {
-                values.put(gaugeName, vs);
-            }
-
+            values.put(gaugeName, vs);
         }
 
         return values;
@@ -176,82 +144,92 @@
      */
     public static double[][] loadDischargeTableValues(
         DischargeTable table,
-        double         scale,
-        boolean        qSorted
+        double         scale
     ) {
         List<DischargeTableValue> dtvs = table.getDischargeTableValues();
 
         final double [][] vs = new double[2][dtvs.size()];
 
-        double lastQ = -Double.MAX_VALUE;
-
         int idx = 0;
         for (DischargeTableValue dtv: dtvs) {
             double q = dtv.getQ().doubleValue();
             vs[0][idx] = q * scale;
             vs[1][idx] = dtv.getW().doubleValue() * scale;
             ++idx;
-
-            if (qSorted && lastQ > q) {
-                qSorted = false;
-            }
-            lastQ = q;
         }
 
         return vs;
     }
 
+    private static final double EPSILON = 1e-5;
 
-    public static double getQForW(double [][] values, double w) {
+    private static final boolean epsEquals(double a, double b) {
+        return Math.abs(a - b) < EPSILON;
+    }
+
+    private static final boolean between(double a, double b, double x) {
+        if (a > b) { double t = a; a = b; b = t; }
+        return x > a && x < b;
+    }
+
+    public static double [] getQsForW(double [][] values, double w) {
 
         boolean debug = log.isDebugEnabled();
 
         if (debug) {
-            log.debug("calculating getQForW(" + w + ")");
-        }
-
-        int index = Arrays.binarySearch(values[1], w);
-        if (index >= 0) {
-            return values[0][index];
-        }
-
-        index = -index - 1; // insert position
-
-        if (index < 1 || index >= values[0].length) {
-            // do not extraploate
-            if (debug) {
-                log.debug("we do not extrapolate: NaN");
-            }
-            return Double.NaN;
+            log.debug("getQsForW: W = " + w);
         }
 
-        double w1 = values[1][index-1];
-        double w2 = values[1][index  ];
-        double q1 = values[0][index-1];
-        double q2 = values[0][index  ];
+        double [] qs = values[0];
+        double [] ws = values[1];
 
-        // q1 = m*w1 + b
-        // q2 = m*w2 + b
-        // q2 - q1 = m*(w2 - w1)
-        // m = (q2 - q1)/(w2 - w1) # w2 != w1
-        // b = q1 - m*w1
+        int N = Math.min(qs.length, ws.length);
 
-        double q;
-        if (w1 == w2) {
-            q = 0.5*(q1 + q2);
+        if (N == 0) {
             if (debug) {
-                log.debug("same w1 and w2: " + w1);
+                log.debug("Q(" + w + ") = []");
+            }
+            return new double [0];
+        }
+
+        TDoubleArrayList outQs = new TDoubleArrayList();
+
+        if (epsEquals(ws[0], w)) {
+            outQs.add(qs[0]);
+        }
+
+        for (int i = 1; i < N; ++i) {
+            if (epsEquals(ws[i], w)) {
+                outQs.add(qs[i]);
+            }
+            else if (between(ws[i-1], ws[i], w)) {
+                double w1 = ws[i-1];
+                double w2 = ws[i];
+                double q1 = qs[i-1];
+                double q2 = qs[i];
+
+                // q1 = m*w1 + b
+                // q2 = m*w2 + b
+                // q2 - q1 = m*(w2 - w1)
+                // m = (q2 - q1)/(w2 - w1) # w2 != w1
+                // b = q1 - m*w1
+                // w1 != w2
+
+                double m = (q2 - q1)/(w2 - w1);
+                double b = q1 - m*w1;
+                double q = w*m + b;
+
+                outQs.add(q);
             }
         }
-        else {
-            double m = (q2 - q1)/(w2 - w1);
-            double b = q1 - m*w1;
-            q = w*m + b;
+
+        double [] result = outQs.toNativeArray();
+
+        if (debug) {
+            log.debug("Q(" + w + ") = " + Arrays.toString(result));
         }
-        if (debug) {
-            log.debug("Q(" + w + ") = " + q);
-        }
-        return q;
+
+        return result;
     }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :