--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Fri Sep 28 12:14:53 2012 +0200
@@ -0,0 +1,1035 @@
+package de.intevation.flys.artifacts.model;
+
+import java.io.Serializable;
+
+import de.intevation.flys.artifacts.math.Linear;
+import de.intevation.flys.artifacts.math.Function;
+
+import java.util.Arrays;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Collections;
+
+import org.apache.log4j.Logger;
+
+import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
+
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+import org.apache.commons.math.ArgumentOutsideDomainException;
+
+import org.apache.commons.math.exception.MathIllegalArgumentException;
+
+import gnu.trove.TDoubleArrayList;
+
+/**
+ * W, Q and km data from database 'wsts' spiced with interpolation algorithms.
+ */
+public class WstValueTable
+implements   Serializable
+{
+    private static Logger log = Logger.getLogger(WstValueTable.class);
+
+    public static final int DEFAULT_Q_STEPS = 500;
+
+    /**
+     * A Column in the table, typically representing one measurement session.
+     */
+    public static final class Column
+    implements                Serializable
+    {
+        protected String name;
+
+        protected QRangeTree qRangeTree;
+
+        public Column() {
+        }
+
+        public Column(String name) {
+            this.name = name;
+        }
+
+        public String getName() {
+            return name;
+        }
+
+        public void setName(String name) {
+            this.name = name;
+        }
+
+    public QRangeTree getQRangeTree() {
+            return qRangeTree;
+        }
+
+        public void setQRangeTree(QRangeTree qRangeTree) {
+            this.qRangeTree = qRangeTree;
+        }
+    } // class Column
+
+    /**
+     * A (weighted) position used for interpolation.
+     */
+    public static final class QPosition {
+
+        protected int    index;
+        protected double weight;
+
+        public QPosition() {
+        }
+
+        public QPosition(int index, double weight) {
+            this.index  = index;
+            this.weight = weight;
+        }
+
+        public QPosition set(int index, double weight) {
+            this.index  = index;
+            this.weight = weight;
+            return this;
+        }
+
+    } // 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.
+     */
+    public static final class Row
+    implements                Serializable, Comparable<Row>
+    {
+        double    km;
+        double [] ws;
+
+        public Row() {
+        }
+
+        public Row(double km) {
+            this.km = km;
+        }
+
+        public Row(double km, double [] ws) {
+            this(km);
+            this.ws = ws;
+        }
+
+        /**
+         * Compare according to place of measurement (km).
+         */
+        public int compareTo(Row other) {
+            double d = km - other.km;
+            if (d < -0.0001) return -1;
+            if (d >  0.0001) return +1;
+            return 0;
+        }
+
+        /**
+         * Interpolate Ws, given Qs and a km.
+         *
+         * @param iqs Given ("input") Qs.
+         * @param ows Resulting ("output") Ws.
+         * @param table Table of which to use data for interpolation.
+         */
+        public void interpolateW(
+            Row           other,
+            double        km,
+            double []     iqs,
+            double []     ows,
+            WstValueTable table,
+            Calculation   errors
+        ) {
+            double kmWeight = Linear.factor(km, this.km, other.km);
+
+            QPosition qPosition = new QPosition();
+
+            for (int i = 0; i < iqs.length; ++i) {
+                if (table.getQPosition(km, iqs[i], qPosition) != null) {
+                    double wt =       getW(qPosition);
+                    double wo = other.getW(qPosition);
+                    if (Double.isNaN(wt) || Double.isNaN(wo)) {
+                        if (errors != null) {
+                            // TODO: I18N
+                            errors.addProblem(
+                                km, "cannot find w for q = " + iqs[i]);
+                        }
+                        ows[i] = Double.NaN;
+                    }
+                    else {
+                        ows[i] = Linear.weight(kmWeight, wt, wo);
+                    }
+                }
+                else {
+                    if (errors != null) {
+                        // TODO: I18N
+                        errors.addProblem(km, "cannot find q = " + iqs[i]);
+                    }
+                    ows[i] = Double.NaN;
+                }
+            }
+        }
+
+
+        public SplineFunction createSpline(
+            WstValueTable table,
+            Calculation   errors
+        ) {
+            int W = ws.length;
+
+            if (W < 1) {
+                if (errors != null) {
+                    // TODO: I18N
+                    errors.addProblem(km, "no ws found");
+                }
+                return null;
+            }
+
+            double [] splineQs = new double[W];
+
+            for (int i = 0; i < W; ++i) {
+                double sq = table.getQIndex(i, km);
+                if (Double.isNaN(sq) && errors != null) {
+                    // TODO: I18N
+                    errors.addProblem(
+                        km, "no q found in " + (i+1) + " column");
+                }
+                splineQs[i] = sq;
+            }
+
+            try {
+                SplineInterpolator interpolator = new SplineInterpolator();
+                PolynomialSplineFunction spline =
+                    interpolator.interpolate(splineQs, ws);
+
+                return new SplineFunction(spline, splineQs, ws);
+            }
+            catch (MathIllegalArgumentException miae) {
+                if (errors != null) {
+                    // TODO: I18N
+                    errors.addProblem(km, "spline creation failed");
+                }
+                log.error("spline creation failed", miae);
+            }
+            return null;
+        }
+
+        public SplineFunction createSpline(
+            Row           other,
+            double        km,
+            WstValueTable table,
+            Calculation   errors
+        ) {
+            int W = Math.min(ws.length, other.ws.length);
+
+            if (W < 1) {
+                if (errors != null) {
+                    // TODO: I18N
+                    errors.addProblem("no ws found");
+                }
+                return null;
+            }
+
+            double factor = Linear.factor(km, this.km, other.km);
+
+            double [] splineQs = new double[W];
+            double [] splineWs = new double[W];
+
+            for (int i = 0; i < W; ++i) {
+                double wws = Linear.weight(factor, ws[i], other.ws[i]);
+                double wqs = Linear.weight(
+                    factor,
+                    table.getQIndex(i, km),
+                    table.getQIndex(i, other.km));
+
+                if (Double.isNaN(wws) || Double.isNaN(wqs)) {
+                    if (errors != null) {
+                        // TODO: I18N
+                        errors.addProblem(km, "cannot find w or q");
+                    }
+                }
+
+                splineWs[i] = wws;
+                splineQs[i] = wqs;
+            }
+
+            SplineInterpolator interpolator = new SplineInterpolator();
+
+            try {
+                PolynomialSplineFunction spline =
+                    interpolator.interpolate(splineQs, splineWs);
+
+                return new SplineFunction(spline, splineQs, splineWs);
+            }
+            catch (MathIllegalArgumentException miae) {
+                if (errors != null) {
+                    // TODO: I18N
+                    errors.addProblem(km, "spline creation failed");
+                }
+                log.error("spline creation failed", miae);
+            }
+
+            return null;
+        }
+
+        public double [][] interpolateWQ(
+            Row           other,
+            double        km,
+            int           steps,
+            WstValueTable table,
+            Calculation   errors
+        ) {
+            SplineFunction sf = createSpline(other, km, table, errors);
+
+            return sf != null
+                ? sf.sample(steps, km, errors)
+                : new double[2][0];
+        }
+
+
+        public double [][] interpolateWQ(
+            int           steps,
+            WstValueTable table,
+            Calculation   errors
+        ) {
+            SplineFunction sf = createSpline(table, errors);
+
+            return sf != null
+                ? sf.sample(steps, km, errors)
+                : new double[2][0];
+        }
+
+
+        public double getW(QPosition qPosition) {
+            int    index  = qPosition.index;
+            double weight = qPosition.weight;
+
+            return weight == 1.0
+                ? ws[index]
+                : Linear.weight(weight, ws[index-1], ws[index]);
+        }
+
+        public double getW(
+            Row       other,
+            double    km,
+            QPosition qPosition
+        ) {
+            double kmWeight = Linear.factor(km, this.km, other.km);
+
+            int    index  = qPosition.index;
+            double weight = qPosition.weight;
+
+            double tw, ow;
+
+            if (weight == 1.0) {
+                tw = ws[index];
+                ow = other.ws[index];
+            }
+            else {
+                tw = Linear.weight(weight, ws[index-1], ws[index]);
+                ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
+            }
+
+            return Linear.weight(kmWeight, tw, ow);
+        }
+
+        public double [] findQsForW(double w, WstValueTable table) {
+
+            TDoubleArrayList qs = new TDoubleArrayList();
+
+            if (ws.length > 0 && Math.abs(ws[0]-w) < 0.000001) {
+                double q = table.getQIndex(0, km);
+                if (!Double.isNaN(q)) {
+                    qs.add(q);
+                }
+            }
+
+            for (int i = 1; i < ws.length; ++i) {
+                double w2 = ws[i];
+                if (Double.isNaN(w2)) {
+                    continue;
+                }
+                if (Math.abs(w2-w) < 0.000001) {
+                    double q = table.getQIndex(i, km);
+                    if (!Double.isNaN(q)) {
+                        qs.add(q);
+                    }
+                    continue;
+                }
+                double w1 = ws[i-1];
+                if (Double.isNaN(w1)) {
+                    continue;
+                }
+
+                if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
+                    continue;
+                }
+
+                double q1 = table.getQIndex(i-1, km);
+                double q2 = table.getQIndex(i,   km);
+                if (Double.isNaN(q1) || Double.isNaN(q2)) {
+                    continue;
+                }
+
+                double q = Linear.linear(w, w1, w2, q1, q2);
+                qs.add(q);
+            }
+
+            return qs.toNativeArray();
+        }
+
+        public double [] findQsForW(
+            Row           other,
+            double        w,
+            double        km, 
+            WstValueTable table
+        ) {
+            TDoubleArrayList qs = new TDoubleArrayList();
+
+            double factor = Linear.factor(km, this.km, other.km);
+
+            if (ws.length > 0) {
+                double wt = Linear.weight(factor, ws[0], other.ws[0]);
+                if (!Double.isNaN(wt)) {
+                    double q = table.getQIndex(0, km);
+                    if (!Double.isNaN(q)) {
+                        qs.add(q);
+                    }
+                }
+            }
+
+            for (int i = 1; i < ws.length; ++i) {
+                double w2 = Linear.weight(factor, ws[i], other.ws[i]);
+                if (Double.isNaN(w2)) {
+                    continue;
+                }
+                if (Math.abs(w2-w) < 0.000001) {
+                    double q = table.getQIndex(i, km);
+                    if (!Double.isNaN(q)) {
+                        qs.add(q);
+                    }
+                    continue;
+                }
+                double w1 = Linear.weight(factor, ws[i-1], other.ws[i-1]);
+                if (Double.isNaN(w1)) {
+                    continue;
+                }
+
+                if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
+                    continue;
+                }
+
+                double q1 = table.getQIndex(i-1, km);
+                double q2 = table.getQIndex(i,   km);
+                if (Double.isNaN(q1) || Double.isNaN(q2)) {
+                    continue;
+                }
+
+                double q = Linear.linear(w, w1, w2, q1, q2);
+                qs.add(q);
+            }
+
+            return qs.toNativeArray();
+        }
+    } // class Row
+
+    /** Rows in table. */
+    protected List<Row> rows;
+
+    /** Columns in table. */
+    protected Column [] columns;
+
+    public WstValueTable() {
+        rows = new ArrayList<Row>();
+    }
+
+    public WstValueTable(Column [] columns) {
+        this();
+        this.columns = columns;
+    }
+
+    public WstValueTable(Column [] columns, List<Row> rows) {
+        this.columns = columns;
+        this.rows    = rows;
+    }
+
+    /** Sort rows (by km). */
+    public void sortRows() {
+        Collections.sort(rows);
+    }
+
+    /**
+     * @param km Given kilometer.
+     * @param qs Given Q values.
+     * @param ws output parameter.
+     */
+    public double [] interpolateW(double km, double [] qs, double [] ws) {
+        return interpolateW(km, qs, ws, null);
+    }
+
+
+    /**
+     * @param ws (output parameter), gets returned.
+     * @return output parameter ws.
+     */
+    public double [] interpolateW(
+        double      km,
+        double []   qs,
+        double []   ws,
+        Calculation errors
+    ) {
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        QPosition qPosition = new QPosition();
+
+        if (rowIndex >= 0) { // direct row match
+            Row row = rows.get(rowIndex);
+            for (int i = 0; i < qs.length; ++i) {
+                if (getQPosition(km, qs[i], qPosition) == null) {
+                    if (errors != null) {
+                        // TODO: I18N
+                        errors.addProblem(km, "cannot find q = " + qs[i]);
+                    }
+                    ws[i] = Double.NaN;
+                }
+                else {
+                    if (Double.isNaN(ws[i] = row.getW(qPosition))
+                    && errors != null) {
+                        // TODO: I18N
+                        errors.addProblem(
+                            km, "cannot find w for q = " + qs[i]);
+                    }
+                }
+            }
+        }
+        else { // needs bilinear interpolation
+            rowIndex = -rowIndex -1;
+
+            if (rowIndex < 1 || rowIndex >= rows.size()) {
+                // do not extrapolate
+                Arrays.fill(ws, Double.NaN);
+                if (errors != null) {
+                    // TODO: I18N
+                    errors.addProblem(km, "km not found");
+                }
+            }
+            else {
+                Row r1 = rows.get(rowIndex-1);
+                Row r2 = rows.get(rowIndex);
+                r1.interpolateW(r2, km, qs, ws, this, errors);
+            }
+        }
+
+        return ws;
+    }
+
+    /**
+     * Interpolate W and Q values at a given km.
+     */
+    public double [][] interpolateWQ(double km) {
+        return interpolateWQ(km, null);
+    }
+
+    /**
+     * Interpolate W and Q values at a given km.
+     */
+    public double [][] interpolateWQ(double km, Calculation errors) {
+        return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
+    }
+
+    /**
+     * Interpolate W and Q values at a given km.
+     */
+    public double [][] interpolateWQ(double km, int steps, Calculation errors) {
+
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        if (rowIndex >= 0) { // direct row match
+            Row row = rows.get(rowIndex);
+            return row.interpolateWQ(steps, this, errors);
+        }
+
+        rowIndex = -rowIndex -1;
+
+        if (rowIndex < 1 || rowIndex >= rows.size()) {
+            // do not extrapolate
+            if (errors != null) {
+                // TODO: I18N
+                errors.addProblem(km, "km not found");
+            }
+            return new double[2][0];
+        }
+
+        Row r1 = rows.get(rowIndex-1);
+        Row r2 = rows.get(rowIndex);
+
+        return r1.interpolateWQ(r2, km, steps, this, errors);
+    }
+
+    public boolean interpolate(
+        double    km,
+        double [] out,
+        QPosition qPosition,
+        Function  qFunction
+    ) {
+        int R1 = rows.size()-1;
+
+        out[1] = qFunction.value(getQ(qPosition, km));
+
+        if (Double.isNaN(out[1])) {
+            return false;
+        }
+
+        QPosition nPosition = new QPosition();
+        if (getQPosition(km, out[1], nPosition) == null) {
+            return false;
+        }
+
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        if (rowIndex >= 0) {
+            // direct row match
+            out[0] = rows.get(rowIndex).getW(nPosition);
+            return !Double.isNaN(out[0]);
+        }
+
+        rowIndex = -rowIndex -1;
+
+        if (rowIndex < 1 || rowIndex > R1) {
+            // do not extrapolate
+            return false;
+        }
+
+        Row r1 = rows.get(rowIndex-1);
+        Row r2 = rows.get(rowIndex);
+        out[0] = r1.getW(r2, km, nPosition);
+
+        return !Double.isNaN(out[0]);
+    }
+
+
+    /**
+     * Look up interpolation of a Q at given positions.
+     *
+     * @param q           the non-interpolated Q value.
+     * @param referenceKm the reference km (e.g. gauge position).
+     * @param kms         positions for which to interpolate.
+     * @param ws          (output) resulting interpolated ws.
+     * @param qs          (output) resulting interpolated qs.
+     * @param errors      calculation object to store errors.
+     */
+    public QPosition interpolate(
+        double      q,
+        double      referenceKm,
+        double []   kms,
+        double []   ws,
+        double []   qs,
+        Calculation errors
+    ) {
+        return interpolate(
+            q, referenceKm, kms, ws, qs, 0, kms.length, errors);
+    }
+
+    public QPosition interpolate(
+        double      q,
+        double      referenceKm,
+        double []   kms,
+        double []   ws,
+        double []   qs,
+        int         startIndex,
+        int         length,
+        Calculation errors
+    ) {
+        QPosition qPosition = getQPosition(referenceKm, q);
+
+        if (qPosition == null) {
+            // we cannot locate q at km
+            Arrays.fill(ws, Double.NaN);
+            Arrays.fill(qs, Double.NaN);
+            if (errors != null) {
+                errors.addProblem(referenceKm, "cannot find q");
+            }
+            return null;
+        }
+
+        Row kmKey = new Row();
+
+        int R1 = rows.size()-1;
+
+        for (int i = startIndex, end = startIndex+length; i < end; ++i) {
+
+            if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
+                if (errors != null) {
+                    errors.addProblem(kms[i], "cannot find q");
+                }
+                ws[i] = Double.NaN;
+                continue;
+            }
+
+            kmKey.km = kms[i];
+            int rowIndex = Collections.binarySearch(rows, kmKey);
+
+            if (rowIndex >= 0) {
+                // direct row match
+                if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
+                && errors != null) {
+                    errors.addProblem(kms[i], "cannot find w");
+                }
+                continue;
+            }
+
+            rowIndex = -rowIndex -1;
+
+            if (rowIndex < 1 || rowIndex > R1) {
+                // do not extrapolate
+                if (errors != null) {
+                    errors.addProblem(kms[i], "cannot find km");
+                }
+                ws[i] = Double.NaN;
+                continue;
+            }
+            Row r1 = rows.get(rowIndex-1);
+            Row r2 = rows.get(rowIndex);
+
+            if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
+            && errors != null) {
+                errors.addProblem(kms[i], "cannot find w");
+            }
+        }
+
+        return qPosition;
+    }
+
+    /**
+     * Linearly interpolate w at a km at a column of two rows.
+     *
+     * @param km   position for which to interpolate.
+     * @param row1 first row.
+     * @param row2 second row.
+     * @param col  column-index at which to look.
+     *
+     * @return Linearly interpolated w, NaN if one of the given rows was null.
+     */
+    public static double linearW(double km, Row row1, Row row2, int col) {
+        if (row1 == null || row2 == null) {
+            return Double.NaN;
+        }
+
+        return Linear.linear(km,
+            row1.km, row2.km,
+            row1.ws[col], row2.ws[col]);
+    }
+
+    /**
+     * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
+     * of other columns).
+     * @param km position (km) at which to interpolate/lookup.
+     * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
+     */
+    public double [][] interpolateWQColumnwise(double km) {
+        log.debug("WstValueTable.interpolateWQColumnwise");
+        double [] qs = new double[columns.length];
+        double [] ws = new double[columns.length];
+
+        // Find out row from where we will start searching.
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        if (rowIndex < 0) {
+            rowIndex = -rowIndex -1;
+        }
+        if (rowIndex >= rows.size()) {
+            rowIndex = rows.size() -1;
+        }
+
+        Row startRow = rows.get(rowIndex);
+
+        for (int col = 0; col < columns.length; col++) {
+            qs[col] = columns[col].getQRangeTree().findQ(km);
+            if (startRow.km == km && startRow.ws[col] != Double.NaN) {
+                // Great. W is defined at km.
+                ws[col] = startRow.ws[col];
+                continue;
+            }
+
+            // Search neighbouring rows that define w at this col.
+            Row rowBefore = null;
+            Row rowAfter  = null;
+            for (int before = rowIndex -1; before >= 0; before--) {
+                if (!Double.isNaN(rows.get(before).ws[col])) {
+                    rowBefore = rows.get(before);
+                    break;
+                }
+            }
+            for (int after = rowIndex, R = rows.size(); after < R; after++) {
+                if (!Double.isNaN(rows.get(after).ws[col])) {
+                    rowAfter = rows.get(after);
+                    break;
+                }
+            }
+
+            ws[col] = linearW(km, rowBefore, rowAfter, col);
+        }
+
+        return new double [][] {qs, ws};
+    }
+
+    public double [] findQsForW(double km, double w) {
+
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
+
+        if (rowIndex >= 0) {
+            return rows.get(rowIndex).findQsForW(w, this);
+        }
+
+        rowIndex = -rowIndex - 1;
+
+        if (rowIndex < 1 || rowIndex >= rows.size()) {
+            // Do not extrapolate.
+            return new double[0];
+        }
+
+        // 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());
+    }
+
+    public QPosition getQPosition(double km, double q, QPosition qPosition) {
+
+        if (columns.length == 0) {
+            return null;
+        }
+
+        double qLast = columns[0].getQRangeTree().findQ(km);
+
+        if (Math.abs(qLast - q) < 0.00001) {
+            return qPosition.set(0, 1d);
+        }
+
+        for (int i = 1; i < columns.length; ++i) {
+            double qCurrent = columns[i].getQRangeTree().findQ(km);
+            if (Math.abs(qCurrent - q) < 0.00001) {
+                return qPosition.set(i, 1d);
+            }
+
+            double qMin, qMax;
+            if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
+            else                  { qMin = qCurrent; qMax = qLast; }
+
+            if (q > qMin && q < qMax) {
+                double weight = Linear.factor(q, qLast, qCurrent);
+                return qPosition.set(i, weight);
+            }
+            qLast = qCurrent;
+        }
+
+        return null;
+    }
+
+    public double getQIndex(int index, double km) {
+        return columns[index].getQRangeTree().findQ(km);
+    }
+
+    public double getQ(QPosition qPosition, double km) {
+        int    index  = qPosition.index;
+        double weight = qPosition.weight;
+
+        if (weight == 1d) {
+            return columns[index].getQRangeTree().findQ(km);
+        }
+        double q1 = columns[index-1].getQRangeTree().findQ(km);
+        double q2 = columns[index  ].getQRangeTree().findQ(km);
+        return Linear.weight(weight, q1, q2);
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :