# HG changeset patch # User Sascha L. Teichmann # Date 1322247131 0 # Node ID 2730d17df0219d03b028d52893f1430198164498 # Parent 1d991c91285be3a912539dec0553c3e15604ebb8 Added the implementation of the 'Bezugslinienverfahren'. flys-artifacts/trunk@3320 c6561f87-3c4e-4783-a992-168aeb5c3f6f diff -r 1d991c91285b -r 2730d17df021 flys-artifacts/ChangeLog --- a/flys-artifacts/ChangeLog Fri Nov 25 17:27:40 2011 +0000 +++ b/flys-artifacts/ChangeLog Fri Nov 25 18:52:11 2011 +0000 @@ -1,3 +1,10 @@ +2011-11-25 Sascha L. Teichmann + + * src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java(relateWs): + Added the implementation of the 'Bezugslinienverfahren'. Should + be complete but needs testing! + TODO: Setup a Calculation and integrate it into WINFO. + 2011-11-25 Sascha L. Teichmann * src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java: diff -r 1d991c91285b -r 2730d17df021 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 Fri Nov 25 17:27:40 2011 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java Fri Nov 25 18:52:11 2011 +0000 @@ -90,6 +90,82 @@ } // class Position + public static final class SplineFunction { + + public PolynomialSplineFunction spline; + public double [] splineQs; + public double [] splineWs; + + public SplineFunction( + PolynomialSplineFunction spline, + double [] splineQs, + double [] splineWs + ) { + this.spline = spline; + this.splineQs = splineQs; + this.splineWs = splineWs; + } + + public double [][] sample( + int numSamples, + double km, + Calculation errors + ) { + double minQ = getQMin(); + double maxQ = getQMax(); + + double [] outWs = new double[numSamples]; + double [] outQs = new double[numSamples]; + + Arrays.fill(outWs, Double.NaN); + Arrays.fill(outQs, Double.NaN); + + double stepWidth = (maxQ - minQ)/numSamples; + + try { + double q = minQ; + for (int i = 0; i < outWs.length; ++i, q += stepWidth) { + outWs[i] = spline.value(outQs[i] = q); + } + } + catch (ArgumentOutsideDomainException aode) { + if (errors != null) { + // TODO: I18N + errors.addProblem(km, "spline interpolation failed"); + } + log.error("spline interpolation failed.", aode); + } + + return new double [][] { outWs, outQs }; + } + + public double getQMin() { + return Math.min(splineQs[0], splineQs[splineQs.length-1]); + } + + public double getQMax() { + return Math.max(splineQs[0], splineQs[splineQs.length-1]); + } + + /** Constructs a continues index between the columns to Qs. */ + public PolynomialSplineFunction createIndexQRelation() { + + double [] indices = new double[splineQs.length]; + for (int i = 0; i < indices.length; ++i) { + indices[i] = i; + } + + try { + SplineInterpolator interpolator = new SplineInterpolator(); + return interpolator.interpolate(indices, splineQs); + } + catch (MathIllegalArgumentException miae) { + // Ignore me! + } + return null; + } + } // class SplineFunction + /** * A row, typically a position where measurements were taken. */ @@ -166,56 +242,6 @@ } } - public static final class SplineFunction { - - public PolynomialSplineFunction spline; - public double [] splineQs; - public double [] splineWs; - - public SplineFunction( - PolynomialSplineFunction spline, - double [] splineQs, - double [] splineWs - ) { - this.spline = spline; - this.splineQs = splineQs; - this.splineWs = splineWs; - } - - public double [][] sample( - int numSamples, - double km, - Calculation errors - ) { - double q1 = splineQs[0], qN = splineQs[splineQs.length-1]; - double minQ = Math.min(q1, qN); - double maxQ = Math.max(q1, qN); - - double [] outWs = new double[numSamples]; - double [] outQs = new double[numSamples]; - - Arrays.fill(outWs, Double.NaN); - Arrays.fill(outQs, Double.NaN); - - double stepWidth = (maxQ - minQ)/numSamples; - - try { - double q = minQ; - for (int i = 0; i < outWs.length; ++i, q += stepWidth) { - outWs[i] = spline.value(outQs[i] = q); - } - } - catch (ArgumentOutsideDomainException aode) { - if (errors != null) { - // TODO: I18N - errors.addProblem(km, "spline interpolation failed"); - } - log.error("spline interpolation failed.", aode); - } - - return new double [][] { outWs, outQs }; - } - } // class SplineFunction public SplineFunction createSpline( WstValueTable table, @@ -821,7 +847,6 @@ public double [] findQsForW(double km, double w) { - int rowIndex = Collections.binarySearch(rows, new Row(km)); if (rowIndex >= 0) { @@ -831,18 +856,130 @@ rowIndex = -rowIndex - 1; if (rowIndex < 1 || rowIndex >= rows.size()) { - // do not extrapolate + // Do not extrapolate. return new double[0]; } - - // Needs bilinear interpolation + // Needs bilinear interpolation. Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); return r1.findQsForW(r2, w, km, this); } + protected SplineFunction createSpline(double km, Calculation errors) { + + int rowIndex = Collections.binarySearch(rows, new Row(km)); + + if (rowIndex >= 0) { + SplineFunction sf = rows.get(rowIndex).createSpline(this, errors); + if (sf == null && errors != null) { + // TODO: I18N + errors.addProblem(km, "cannot create w/q relation"); + } + return sf; + } + + rowIndex = -rowIndex - 1; + + if (rowIndex < 1 || rowIndex >= rows.size()) { + // Do not extrapolate. + if (errors != null) { + // TODO: I18N + errors.addProblem(km, "km not found"); + } + return null; + } + + // Needs bilinear interpolation. + Row r1 = rows.get(rowIndex-1); + Row r2 = rows.get(rowIndex); + + SplineFunction sf = r1.createSpline(r2, km, this, errors); + if (sf == null && errors != null) { + // TODO: I18N + errors.addProblem(km, "cannot create w/q relation"); + } + + return sf; + } + + /** 'Bezugslinienverfahren' */ + public double [][] relateWs( + double km1, + double km2, + int numSamples, + Calculation errors + ) { + SplineFunction sf1 = createSpline(km1, errors); + if (sf1 == null) { + return new double[2][0]; + } + + SplineFunction sf2 = createSpline(km1, errors); + if (sf2 == null) { + return new double[2][0]; + } + + PolynomialSplineFunction iQ1 = sf1.createIndexQRelation(); + if (iQ1 == null) { + if (errors != null) { + // TODO: I18N + errors.addProblem(km1, "cannot create index/q relation"); + } + return new double[2][0]; + } + + PolynomialSplineFunction iQ2 = sf1.createIndexQRelation(); + if (iQ2 == null) { + if (errors != null) { + // TODO: I18N + errors.addProblem(km2, "cannot create index/q relation"); + } + return new double[2][0]; + } + + int N = sf1.splineQs.length; + double stepWidth = N/(double)numSamples; + + PolynomialSplineFunction qW1 = sf1.spline; + PolynomialSplineFunction qW2 = sf2.spline; + + double [] ws1 = new double[numSamples]; + double [] ws2 = new double[numSamples]; + + Arrays.fill(ws1, Double.NaN); + Arrays.fill(ws2, Double.NaN); + + boolean hadErrors = false; + + double p = 0d; + for (int i = 0; i < numSamples; ++i, p += stepWidth) { + try { + double q1 = iQ1.value(p); + double w1 = qW1.value(q1); + double q2 = iQ2.value(p); + double w2 = qW2.value(q2); + ws1[i] = w1; + ws2[i] = w2; + } + catch (ArgumentOutsideDomainException aode) { + if (!hadErrors) { + // XXX: I'm not sure if this really can happen + // and if we should report this more than once. + hadErrors = true; + if (errors != null) { + // TODO: I18N + log.debug("W~W failed", aode); + errors.addProblem("W~W failed"); + } + } + } + } + + return new double [][] { ws1, ws2 }; + } + public QPosition getQPosition(double km, double q) { return getQPosition(km, q, new QPosition()); }