# HG changeset patch # User Sascha L. Teichmann # Date 1303235019 0 # Node ID 64cfbd631f294041e9fc8142ff61fc2a4094467b # Parent b7c8df643dc45741d87094c7b4a59b36cdfab028 Implemented the "Wasserstand/Wasserspiegellagen" calculation. flys-artifacts/trunk@1733 c6561f87-3c4e-4783-a992-168aeb5c3f6f diff -r b7c8df643dc4 -r 64cfbd631f29 flys-artifacts/ChangeLog --- 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 + + * 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 * src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java: diff -r b7c8df643dc4 -r 64cfbd631f29 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.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 { 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