changeset 1939:2730d17df021

Added the implementation of the 'Bezugslinienverfahren'. flys-artifacts/trunk@3320 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 25 Nov 2011 18:52:11 +0000
parents 1d991c91285b
children 0d12e70766c8
files flys-artifacts/ChangeLog flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java
diffstat 2 files changed, 198 insertions(+), 54 deletions(-) [+]
line wrap: on
line diff
--- 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	<sascha.teichmann@intevation.de>
+
+	* 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	<sascha.teichmann@intevation.de>
 
 	* 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());
     }

http://dive4elements.wald.intevation.org