changeset 335:64cfbd631f29

Implemented the "Wasserstand/Wasserspiegellagen" calculation. flys-artifacts/trunk@1733 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 19 Apr 2011 17:43:39 +0000
parents b7c8df643dc4
children 7f13ed751277
files flys-artifacts/ChangeLog flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java
diffstat 2 files changed, 187 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/flys-artifacts/ChangeLog	Tue Apr 19 17:37:21 2011 +0000
+++ b/flys-artifacts/ChangeLog	Tue Apr 19 17:43:39 2011 +0000
@@ -1,3 +1,16 @@
+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
+	  q values and returns an equal sized array of w values.
+	  This is essentially the "Wasserstand/Wasserspiegellagen" calculation
+	  of desktop FLYS.
+
+	  If you want to do a calculation with given w values you have
+	  to convert the w values with DischargeTables.getQForW() first.
+
+	  !!! This code needs heavy testing !!!
+
 2011-04-19	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java:
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Tue Apr 19 17:37:21 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Tue Apr 19 17:43:39 2011 +0000
@@ -8,6 +8,7 @@
 
 import de.intevation.flys.backend.SessionHolder;
 
+import java.util.Arrays;
 import java.util.ArrayList;
 import java.util.Comparator;
 import java.util.List;
@@ -53,17 +54,139 @@
     // class Column
 
     public static class Row
-    implements          Serializable
+    implements          Serializable, Comparable<Row>
     {
         double    km;
-        double [] wq;
+        double [] ws;
+        double [] qs;
+        boolean   qSorted;
 
         public Row() {
         }
 
-        public Row(double km, double [] wq) {
+        public Row(double km) {
             this.km = km;
-            this.wq = wq;
+        }
+
+        public Row(double km, double [] ws, double [] qs) {
+            this(km);
+            this.ws = ws;
+            this.qs = qs;
+        }
+
+        public static final double linear(
+            double x,
+            double x1, double x2,
+            double y1, double y2
+        ) {
+            // y1 = m*x1 + b
+            // y2 = m*x2 + b
+            // y2 - y1 = m*(x2 - x1)
+            // m = (y2 - y1)/(x2 - x1) # x2 != x1
+            // b = y1 - m*x1
+
+            if (x1 == x2) {
+                return 0.5*(y1 + y2);
+            }
+            double m = (y2 - y1)/(x2 - x1);
+            double b = y1 - m*x1;
+            return x*m + b;
+        }
+
+        public static final double factor(double x, double p1, double p2) {
+            // 0 = m*p1 + b <=> b = -m*p1
+            // 1 = m*p2 + b
+            // 1 = m*(p2 - p1)
+            // m = 1/(p2 - p1) # p1 != p2
+            // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
+
+            return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
+        }
+
+        public static final double weight(double factor, double a, double b) {
+            return (1.0-factor)*a + factor*b;
+        }
+
+        public double getWForKM(Row other, int index, double km) {
+            double w1  = ws[index];
+            double w2  = other.ws[index];
+            double km1 = this.km;
+            double km2 = other.km;
+            return linear(km, km1, km2, w1, w2);
+        }
+
+        public void interpolateW(
+            Row       other,
+            double    km,
+            double [] input,
+            double [] output
+        ) {
+            double factor = factor(km, this.km, other.km);
+
+            for (int i = 0; i < input.length; ++i) {
+                double q = input[i];
+                int idx1 =       getQIndex(q);
+                int idx2 = other.getQIndex(q);
+
+                double w1 = idx1 >= 0
+                    ? ws[idx1]
+                    : interpolateW(-idx1-1, q);
+
+                double w2 = idx2 >= 0
+                    ? other.ws[idx2]
+                    : other.interpolateW(-idx2-1, q);
+
+                output[i] = weight(factor, w1, w2);
+            }
+        }
+
+        public double interpolateW(int idx, double q) {
+            return idx < 1 || idx >= qs.length
+                ? Double.NaN // do not extrapolate
+                : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]);
+        }
+
+        public double interpolateW(double q) {
+            if (Double.isNaN(q)) return Double.NaN;
+            int index = getQIndex(q);
+            return index >= 0 ? ws[index] : interpolateW(-index -1);
+        }
+
+        public int getQIndex(double q) {
+            return qSorted ? binaryQIndex(q) : linearQIndex(q);
+        }
+
+        public int binaryQIndex(double q) {
+            return Arrays.binarySearch(qs, q);
+        }
+
+        public int linearQIndex(double q) {
+            switch (qs.length) {
+                case 0: break;
+                case 1:
+                    if (qs[0] == q) return 0;
+                    if (qs[0] <  q) return -(1+1);
+                    return -(0+1);
+                default:
+                    for (int i = 1; i < qs.length; ++i) {
+                        double qa = qs[i-1];
+                        double qb = qs[i];
+                        if (qa == q) return i-1;
+                        if (qb == q) return i;
+                        if (qa > qb) { double t = qa; qa = qb; qb = t; }
+                        if (q > qa && q < qb) return -(i+1);
+                    }
+                    return -qs.length;
+            }
+
+            return -1;
+        }
+
+        public int compareTo(Row other) {
+            double d = km - other.km;
+            if (d < 0.0) return -1;
+            if (d > 0.0) return +1;
+            return 0;
         }
     }
     // class Row
@@ -81,6 +204,36 @@
         this.columns = columns;
     }
 
+
+    public double [] interpolateW(double km, double [] qs) {
+
+        double [] ws = new double[qs.length];
+
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        if (rowIndex >= 0) { // direct row match
+            Row row = rows.get(rowIndex);
+            for (int i = 0; i < qs.length; ++i) {
+                ws[i] = row.interpolateW(qs[i]);
+            }
+        }
+        else { // needs bilinear interpolation
+            rowIndex = -rowIndex -1;
+
+            if (rowIndex < 1 || rowIndex >= rows.size()) {
+                // do not extrapolate
+                Arrays.fill(ws, Double.NaN);
+            }
+            else {
+                rows.get(rowIndex-1).interpolateW(
+                    rows.get(rowIndex),
+                    km, qs, ws);
+            }
+        }
+
+        return ws;
+    }
+
     public static WstValueTable getTable(River river) {
         return getTable(river, 0);
     }
@@ -133,6 +286,9 @@
         int lastColumnNo = -1;
         Row row          = null;
 
+        Double lastQ     = -Double.MAX_VALUE;
+        boolean qSorted  = true;
+
         for (Iterator<Object []> iter = sqlQuery.iterate(); iter.hasNext();) {
             Object [] result = iter.next();
             double km    = (Double) result[0];
@@ -142,14 +298,24 @@
 
             if (columnNo > lastColumnNo) { // new row
                 if (row != null) {
+                    row.qSorted = qSorted;
                     valueTable.rows.add(row);
                 }
-                row = new Row(km, new double[(columnNo+1) << 1]);
+                row = new Row(
+                    km,
+                    new double[columnNo+1],
+                    new double[columnNo+1]);
+                lastQ = -Double.MAX_VALUE;
+                qSorted = true;
             }
-            int idx = columnNo << 1;
 
-            row.wq[idx  ] = w != null ? w : Double.NaN;
-            row.wq[idx+1] = q != null ? q : Double.NaN;
+            row.ws[columnNo] = w != null ? w : Double.NaN;
+            row.qs[columnNo] = q != null ? q : Double.NaN;
+
+            if (qSorted && (q == null || lastQ > q)) {
+                qSorted = false;
+            }
+            lastQ = q;
 
             lastColumnNo = columnNo;
         }

http://dive4elements.wald.intevation.org