# HG changeset patch # User Sascha L. Teichmann # Date 1305493721 0 # Node ID 909196be11a0622e9207b7b8884865145acd7228 # Parent c6b7ca0febcb1665896d8fdd467649afc5b4988b First step to calculate "W fuer ungleichwertige Abfluesse" correctly. flys/issue55 flys-artifacts/trunk@1925 c6561f87-3c4e-4783-a992-168aeb5c3f6f diff -r c6b7ca0febcb -r 909196be11a0 flys-artifacts/ChangeLog --- a/flys-artifacts/ChangeLog Sun May 15 10:59:08 2011 +0000 +++ b/flys-artifacts/ChangeLog Sun May 15 21:08:41 2011 +0000 @@ -1,3 +1,16 @@ +2011-05-10 Sascha L. Teichmann + + First step to calculate "W fuer ungleichwertige Abfluesse" correctly. + flys/issue55 + + * src/main/java/de/intevation/flys/artifacts/math/LinearRemap.java: + New. Remaps "gleichwertige" Q values to the corresponding + "ungleichwertige" Q values depending on km. + + * src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java: + Remap the Q values "ungleichwertig" depending on the + "gleichwertige" ones. + 2011-05-10 Sascha L. Teichmann First step to fix flys/issue69 diff -r c6b7ca0febcb -r 909196be11a0 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/LinearRemap.java --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/LinearRemap.java Sun May 15 21:08:41 2011 +0000 @@ -0,0 +1,78 @@ +package de.intevation.flys.artifacts.math; + +public class LinearRemap +{ + public static class Segment { + + protected Segment next; + + protected double from; + protected double to; + + protected double m; + protected double b; + + public Segment() { + } + + public Segment( + double from, double to, + double m, double b, + Segment next + ) { + this.from = from; + this.to = to; + this.m = m; + this.b = b; + } + + public double eval(double x) { + return m*x + b; + } + } // class Segment + + protected Segment head; + + public LinearRemap() { + } + + public void add( + double from, double to, + double x1, double y1, + double x2, double y2 + ) { + // y1 = m*x1 + b <=> b = y1 - m*x1 + // y2 = m*x2 + b + // y2 - y1 = m*(x2 - x1) + // m = (y2 - y1)/(x2 - x1) + + double m, b; + + if (x2 == x1) { + m = 0.0; + b = 0.5*(y2 + y1); + } + else { + m = (y2 - y1)/(x2 - x1); + b = y1 - m*x1; + } + + head = new Segment(from, to, m, b, head); + } + + public double eval(double pos, double x) { + Segment current = head; + + while (current != null) { + + if (pos >= current.from && pos <= current.to) { + return current.eval(x); + } + + current = current.next; + } + + return Double.NaN; + } +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : diff -r c6b7ca0febcb -r 909196be11a0 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 Sun May 15 10:59:08 2011 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java Sun May 15 21:08:41 2011 +0000 @@ -6,6 +6,8 @@ import de.intevation.flys.model.Wst; import de.intevation.flys.model.WstColumn; +import de.intevation.flys.artifacts.math.LinearRemap; + import de.intevation.flys.artifacts.cache.CacheFactory; import de.intevation.flys.backend.SessionHolder; @@ -441,6 +443,60 @@ oqs[index] = weight(kmWeight, oq, tq); ows[index] = weight(kmWeight, ow, tw); } + + public void storeWQ( + QPosition qPosition, + LinearRemap remap, + double [] ows, + double [] oqs, + int index + ) { + int qIdx = qPosition.index; + double qWeight = qPosition.weight; + + oqs[index] = remap.eval(km, qWeight == 1.0 + ? qs[qIdx] + : weight(qWeight, qs[qIdx-1], qs[qIdx])); + + ows[index] = interpolateW(oqs[index]); + } + + public void storeWQ( + Row other, + double km, + QPosition qPosition, + LinearRemap remap, + double [] ows, + double [] oqs, + int index + ) { + double kmWeight = factor(km, this.km, other.km); + + double tq, tw; + double oq, ow; + + int qIdx = qPosition.index; + double qWeight = qPosition.weight; + + if (qWeight == 1.0) { + tq = remap.eval(this.km, qs[qIdx]); + oq = remap.eval(other.km, other.qs[qIdx]); + } + else { + tq = remap.eval( + this.km, + weight(qWeight, qs[qIdx-1], qs[qIdx])); + oq = remap.eval( + other.km, + weight(qWeight, other.qs[qIdx-1], other.qs[qIdx])); + } + + tw = interpolateW(tq); + ow = other.interpolateW(oq); + + oqs[index] = weight(kmWeight, oq, tq); + ows[index] = weight(kmWeight, ow, tw); + } } // class Row @@ -515,7 +571,40 @@ return r1.interpolateWQ(r2, km, steps); } - public boolean interpolate( + public void interpolate( + double [] kms, + double [] ws, + double [] qs, + QPosition qPosition, + LinearRemap remap + ) { + int R1 = rows.size()-1; + + Row kmKey = new Row(); + for (int i = 0; i < kms.length; ++i) { + kmKey.km = kms[i]; + int rowIndex = Collections.binarySearch(rows, kmKey); + + if (rowIndex >= 0) { + // direct row match + rows.get(rowIndex).storeWQ(qPosition, remap, ws, qs, i); + continue; + } + + rowIndex = -rowIndex -1; + + if (rowIndex < 1 || rowIndex >= R1) { + // do not extrapolate + ws[i] = qs[i] = Double.NaN; + continue; + } + Row r1 = rows.get(rowIndex-1); + Row r2 = rows.get(rowIndex); + r1.storeWQ(r2, kms[i], qPosition, remap, ws, qs, i); + } + } + + public QPosition interpolate( double q, int referenceIndex, double [] kms, @@ -538,7 +627,7 @@ if (rowIndex < 1 || rowIndex >= R1) { // do not extrapolate - return false; + return null; } // interpolation case @@ -550,7 +639,7 @@ if (qPosition == null) { // we cannot locate q in reference row - return false; + return null; } for (int i = 0; i < kms.length; ++i) { @@ -575,7 +664,7 @@ r1.storeWQ(r2, kms[i], qPosition, ws, qs, i); } - return true; + return qPosition; } public static WstValueTable getTable(River river) {