sascha@326: package de.intevation.flys.artifacts.model;
sascha@326: 
sascha@326: import java.io.Serializable;
sascha@326: 
sascha@655: import de.intevation.flys.artifacts.math.Linear;
sascha@655: import de.intevation.flys.artifacts.math.Function;
sascha@427: 
sascha@335: import java.util.Arrays;
sascha@326: import java.util.ArrayList;
sascha@326: import java.util.List;
sascha@326: import java.util.Collections;
sascha@326: 
sascha@339: import org.apache.log4j.Logger;
sascha@339: 
sascha@339: import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
sascha@339: 
sascha@339: import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
sascha@339: 
sascha@339: import org.apache.commons.math.ArgumentOutsideDomainException;
sascha@339: 
ingo@686: import org.apache.commons.math.exception.MathIllegalArgumentException;
ingo@686: 
sascha@1937: import gnu.trove.TDoubleArrayList;
sascha@1937: 
felix@1838: /**
felix@1838:  * W, Q and km data from database 'wsts' spiced with interpolation algorithms.
felix@1838:  */
sascha@326: public class WstValueTable
sascha@326: implements   Serializable
sascha@326: {
sascha@339:     private static Logger log = Logger.getLogger(WstValueTable.class);
sascha@326: 
sascha@395:     public static final int DEFAULT_Q_STEPS = 500;
sascha@339: 
sascha@2253:     public static final int RELATE_WS_SAMPLES = 200;
sascha@2186: 
felix@1721:     /**
felix@1721:      * A Column in the table, typically representing one measurement session.
felix@1721:      */
sascha@633:     public static final class Column
sascha@633:     implements                Serializable
sascha@326:     {
sascha@326:         protected String name;
sascha@326: 
sascha@632:         protected QRangeTree qRangeTree;
sascha@632: 
sascha@326:         public Column() {
sascha@326:         }
sascha@326: 
sascha@326:         public Column(String name) {
sascha@326:             this.name = name;
sascha@326:         }
sascha@326: 
sascha@326:         public String getName() {
sascha@326:             return name;
sascha@326:         }
sascha@326: 
sascha@326:         public void setName(String name) {
sascha@326:             this.name = name;
sascha@326:         }
sascha@632: 
sascha@1907:     public QRangeTree getQRangeTree() {
sascha@632:             return qRangeTree;
sascha@632:         }
sascha@632: 
sascha@632:         public void setQRangeTree(QRangeTree qRangeTree) {
sascha@632:             this.qRangeTree = qRangeTree;
sascha@632:         }
sascha@633:     } // class Column
sascha@326: 
felix@1721:     /**
felix@1721:      * A (weighted) position used for interpolation.
felix@1721:      */
sascha@633:     public static final class QPosition {
sascha@426: 
sascha@426:         protected int    index;
sascha@426:         protected double weight;
sascha@426: 
sascha@426:         public QPosition() {
sascha@426:         }
sascha@426: 
sascha@633:         public QPosition(int index, double weight) {
sascha@426:             this.index  = index;
sascha@426:             this.weight = weight;
sascha@426:         }
sascha@426: 
sascha@633:         public QPosition set(int index, double weight) {
sascha@633:             this.index  = index;
sascha@633:             this.weight = weight;
sascha@633:             return this;
sascha@633:         }
sascha@633: 
sascha@426:     } // class Position
sascha@426: 
sascha@1939:     public static final class SplineFunction {
sascha@1939: 
sascha@1939:         public PolynomialSplineFunction spline;
sascha@1939:         public double []                splineQs;
sascha@1939:         public double []                splineWs;
sascha@1939: 
sascha@1939:         public SplineFunction(
sascha@1939:             PolynomialSplineFunction spline,
sascha@3076:             double []                splineQs,
sascha@1939:             double []                splineWs
sascha@1939:         ) {
sascha@1939:             this.spline   = spline;
sascha@1939:             this.splineQs = splineQs;
sascha@1939:             this.splineWs = splineWs;
sascha@1939:         }
sascha@1939: 
sascha@1939:         public double [][] sample(
sascha@3076:             int         numSamples,
sascha@3076:             double      km,
sascha@1939:             Calculation errors
sascha@1939:         ) {
sascha@1939:             double minQ = getQMin();
sascha@1939:             double maxQ = getQMax();
sascha@1939: 
sascha@1939:             double [] outWs = new double[numSamples];
sascha@1939:             double [] outQs = new double[numSamples];
sascha@1939: 
sascha@1939:             Arrays.fill(outWs, Double.NaN);
sascha@1939:             Arrays.fill(outQs, Double.NaN);
sascha@1939: 
sascha@1939:             double stepWidth = (maxQ - minQ)/numSamples;
sascha@1939: 
sascha@1939:             try {
sascha@1939:                 double q = minQ;
sascha@1939:                 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
sascha@1939:                     outWs[i] = spline.value(outQs[i] = q);
sascha@1939:                 }
sascha@1939:             }
sascha@1939:             catch (ArgumentOutsideDomainException aode) {
sascha@1939:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "spline.interpolation.failed");
sascha@1939:                 }
sascha@1939:                 log.error("spline interpolation failed.", aode);
sascha@1939:             }
sascha@1939: 
sascha@1939:             return new double [][] { outWs, outQs };
sascha@1939:         }
sascha@1939: 
sascha@1939:         public double getQMin() {
sascha@1939:             return Math.min(splineQs[0], splineQs[splineQs.length-1]);
sascha@1939:         }
sascha@1939: 
sascha@1939:         public double getQMax() {
sascha@1939:             return Math.max(splineQs[0], splineQs[splineQs.length-1]);
sascha@1939:         }
sascha@1939: 
sascha@1939:         /** Constructs a continues index between the columns to Qs. */
sascha@1939:         public PolynomialSplineFunction createIndexQRelation() {
sascha@1939: 
sascha@1939:             double [] indices = new double[splineQs.length];
sascha@1939:             for (int i = 0; i < indices.length; ++i) {
sascha@1939:                 indices[i] = i;
sascha@1939:             }
sascha@1939: 
sascha@1939:             try {
sascha@1939:                 SplineInterpolator interpolator = new SplineInterpolator();
sascha@1939:                 return interpolator.interpolate(indices, splineQs);
sascha@1939:             }
sascha@1939:             catch (MathIllegalArgumentException miae) {
sascha@1939:                 // Ignore me!
sascha@1939:             }
sascha@1939:             return null;
sascha@1939:         }
sascha@1939:     } // class SplineFunction
sascha@1939: 
felix@1721:     /**
felix@1838:      * A row, typically a position where measurements were taken.
felix@1721:      */
sascha@633:     public static final class Row
sascha@633:     implements                Serializable, Comparable<Row>
sascha@326:     {
sascha@326:         double    km;
sascha@335:         double [] ws;
sascha@326: 
sascha@326:         public Row() {
sascha@326:         }
sascha@326: 
sascha@335:         public Row(double km) {
sascha@326:             this.km = km;
sascha@335:         }
sascha@335: 
sascha@633:         public Row(double km, double [] ws) {
sascha@335:             this(km);
sascha@335:             this.ws = ws;
sascha@335:         }
sascha@335: 
felix@1838:         /**
felix@1838:          * Compare according to place of measurement (km).
felix@1838:          */
sascha@335:         public int compareTo(Row other) {
sascha@335:             double d = km - other.km;
sascha@458:             if (d < -0.0001) return -1;
sascha@458:             if (d >  0.0001) return +1;
sascha@335:             return 0;
sascha@326:         }
sascha@339: 
felix@1887:         /**
felix@1887:          * Interpolate Ws, given Qs and a km.
felix@1887:          *
felix@1887:          * @param iqs Given ("input") Qs.
felix@1887:          * @param ows Resulting ("output") Ws.
felix@1887:          * @param table Table of which to use data for interpolation.
felix@1887:          */
sascha@633:         public void interpolateW(
sascha@633:             Row           other,
sascha@633:             double        km,
sascha@633:             double []     iqs,
sascha@633:             double []     ows,
ingo@686:             WstValueTable table,
ingo@686:             Calculation   errors
sascha@633:         ) {
sascha@655:             double kmWeight = Linear.factor(km, this.km, other.km);
sascha@339: 
sascha@633:             QPosition qPosition = new QPosition();
sascha@633: 
sascha@633:             for (int i = 0; i < iqs.length; ++i) {
sascha@633:                 if (table.getQPosition(km, iqs[i], qPosition) != null) {
sascha@633:                     double wt =       getW(qPosition);
sascha@633:                     double wo = other.getW(qPosition);
ingo@686:                     if (Double.isNaN(wt) || Double.isNaN(wo)) {
ingo@686:                         if (errors != null) {
ingo@686:                             errors.addProblem(
sascha@2166:                                 km, "cannot.find.w.for.q", iqs[i]);
ingo@686:                         }
ingo@686:                         ows[i] = Double.NaN;
ingo@686:                     }
ingo@686:                     else {
ingo@686:                         ows[i] = Linear.weight(kmWeight, wt, wo);
ingo@686:                     }
sascha@633:                 }
sascha@633:                 else {
ingo@686:                     if (errors != null) {
sascha@2166:                         errors.addProblem(km, "cannot.find.q", iqs[i]);
ingo@686:                     }
sascha@633:                     ows[i] = Double.NaN;
sascha@633:                 }
sascha@633:             }
sascha@633:         }
sascha@633: 
sascha@1938: 
sascha@1938:         public SplineFunction createSpline(
sascha@1938:             WstValueTable table,
sascha@1938:             Calculation   errors
sascha@1938:         ) {
sascha@1938:             int W = ws.length;
sascha@1938: 
sascha@1938:             if (W < 1) {
sascha@1938:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "no.ws.found");
sascha@1938:                 }
sascha@1938:                 return null;
sascha@1938:             }
sascha@1938: 
sascha@1938:             double [] splineQs = new double[W];
sascha@1938: 
sascha@1938:             for (int i = 0; i < W; ++i) {
sascha@1938:                 double sq = table.getQIndex(i, km);
sascha@1938:                 if (Double.isNaN(sq) && errors != null) {
sascha@1938:                     errors.addProblem(
sascha@2166:                         km, "no.q.found.in.column", (i+1));
sascha@1938:                 }
sascha@1938:                 splineQs[i] = sq;
sascha@1938:             }
sascha@1938: 
sascha@1938:             try {
sascha@1938:                 SplineInterpolator interpolator = new SplineInterpolator();
sascha@1938:                 PolynomialSplineFunction spline =
sascha@1938:                     interpolator.interpolate(splineQs, ws);
sascha@1938: 
sascha@1938:                 return new SplineFunction(spline, splineQs, ws);
sascha@1938:             }
sascha@1938:             catch (MathIllegalArgumentException miae) {
sascha@1938:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "spline.creation.failed");
sascha@1938:                 }
sascha@1938:                 log.error("spline creation failed", miae);
sascha@1938:             }
sascha@1938:             return null;
sascha@1938:         }
sascha@1938: 
sascha@1938:         public SplineFunction createSpline(
sascha@633:             Row           other,
sascha@633:             double        km,
ingo@686:             WstValueTable table,
ingo@686:             Calculation   errors
sascha@633:         ) {
sascha@395:             int W = Math.min(ws.length, other.ws.length);
sascha@339: 
sascha@339:             if (W < 1) {
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem("no.ws.found");
ingo@686:                 }
sascha@1938:                 return null;
sascha@339:             }
sascha@339: 
sascha@655:             double factor = Linear.factor(km, this.km, other.km);
sascha@339: 
sascha@1938:             double [] splineQs = new double[W];
sascha@1938:             double [] splineWs = new double[W];
sascha@339: 
sascha@339:             for (int i = 0; i < W; ++i) {
sascha@655:                 double wws = Linear.weight(factor, ws[i], other.ws[i]);
sascha@655:                 double wqs = Linear.weight(
sascha@633:                     factor,
sascha@633:                     table.getQIndex(i, km),
sascha@633:                     table.getQIndex(i, other.km));
sascha@633: 
ingo@686:                 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
ingo@686:                     if (errors != null) {
sascha@2166:                         errors.addProblem(km, "cannot.find.w.or.q");
ingo@686:                     }
ingo@686:                 }
ingo@686: 
sascha@1938:                 splineWs[i] = wws;
sascha@1938:                 splineQs[i] = wqs;
sascha@339:             }
sascha@339: 
sascha@339:             SplineInterpolator interpolator = new SplineInterpolator();
ingo@686: 
ingo@686:             try {
sascha@1938:                 PolynomialSplineFunction spline =
sascha@1938:                     interpolator.interpolate(splineQs, splineWs);
sascha@1938: 
sascha@1938:                 return new SplineFunction(spline, splineQs, splineWs);
ingo@686:             }
ingo@686:             catch (MathIllegalArgumentException miae) {
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "spline.creation.failed");
ingo@686:                 }
felix@1860:                 log.error("spline creation failed", miae);
ingo@686:             }
sascha@339: 
sascha@1938:             return null;
sascha@1938:         }
ingo@686: 
sascha@1938:         public double [][] interpolateWQ(
sascha@1938:             Row           other,
sascha@1938:             double        km,
sascha@1938:             int           steps,
sascha@1938:             WstValueTable table,
sascha@1938:             Calculation   errors
sascha@1938:         ) {
sascha@1938:             SplineFunction sf = createSpline(other, km, table, errors);
sascha@339: 
sascha@1938:             return sf != null
sascha@1938:                 ? sf.sample(steps, km, errors)
sascha@1938:                 : new double[2][0];
sascha@339:         }
sascha@339: 
felix@1721: 
ingo@686:         public double [][] interpolateWQ(
sascha@742:             int           steps,
ingo@686:             WstValueTable table,
ingo@686:             Calculation   errors
ingo@686:         ) {
sascha@1938:             SplineFunction sf = createSpline(table, errors);
sascha@339: 
sascha@1938:             return sf != null
sascha@1938:                 ? sf.sample(steps, km, errors)
sascha@1938:                 : new double[2][0];
sascha@339:         }
sascha@339: 
sascha@426: 
sascha@633:         public double getW(QPosition qPosition) {
sascha@633:             int    index  = qPosition.index;
sascha@633:             double weight = qPosition.weight;
sascha@426: 
sascha@633:             return weight == 1.0
sascha@633:                 ? ws[index]
sascha@655:                 : Linear.weight(weight, ws[index-1], ws[index]);
sascha@426:         }
sascha@426: 
sascha@633:         public double getW(
sascha@426:             Row       other,
sascha@426:             double    km,
sascha@633:             QPosition qPosition
sascha@426:         ) {
sascha@655:             double kmWeight = Linear.factor(km, this.km, other.km);
sascha@426: 
sascha@633:             int    index  = qPosition.index;
sascha@633:             double weight = qPosition.weight;
sascha@426: 
sascha@633:             double tw, ow;
sascha@426: 
sascha@633:             if (weight == 1.0) {
sascha@633:                 tw = ws[index];
sascha@633:                 ow = other.ws[index];
sascha@426:             }
sascha@426:             else {
sascha@742:                 tw = Linear.weight(weight, ws[index-1], ws[index]);
sascha@742:                 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
sascha@426:             }
sascha@426: 
sascha@655:             return Linear.weight(kmWeight, tw, ow);
sascha@427:         }
sascha@1937: 
sascha@1937:         public double [] findQsForW(double w, WstValueTable table) {
sascha@1937: 
sascha@1937:             TDoubleArrayList qs = new TDoubleArrayList();
sascha@1937: 
sascha@1937:             if (ws.length > 0 && Math.abs(ws[0]-w) < 0.000001) {
sascha@1937:                 double q = table.getQIndex(0, km);
sascha@1937:                 if (!Double.isNaN(q)) {
sascha@1937:                     qs.add(q);
sascha@1937:                 }
sascha@1937:             }
sascha@1937: 
sascha@1937:             for (int i = 1; i < ws.length; ++i) {
sascha@1937:                 double w2 = ws[i];
sascha@1937:                 if (Double.isNaN(w2)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937:                 if (Math.abs(w2-w) < 0.000001) {
sascha@1937:                     double q = table.getQIndex(i, km);
sascha@1937:                     if (!Double.isNaN(q)) {
sascha@1937:                         qs.add(q);
sascha@1937:                     }
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937:                 double w1 = ws[i-1];
sascha@1937:                 if (Double.isNaN(w1)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937: 
sascha@1937:                 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937: 
sascha@1937:                 double q1 = table.getQIndex(i-1, km);
sascha@1937:                 double q2 = table.getQIndex(i,   km);
sascha@1937:                 if (Double.isNaN(q1) || Double.isNaN(q2)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937: 
sascha@1937:                 double q = Linear.linear(w, w1, w2, q1, q2);
sascha@1937:                 qs.add(q);
sascha@1937:             }
sascha@1937: 
sascha@1937:             return qs.toNativeArray();
sascha@1937:         }
sascha@1937: 
sascha@1937:         public double [] findQsForW(
sascha@1937:             Row           other,
sascha@1937:             double        w,
sascha@3076:             double        km,
sascha@1937:             WstValueTable table
sascha@1937:         ) {
sascha@1937:             TDoubleArrayList qs = new TDoubleArrayList();
sascha@1937: 
sascha@1937:             double factor = Linear.factor(km, this.km, other.km);
sascha@1937: 
sascha@1937:             if (ws.length > 0) {
sascha@1937:                 double wt = Linear.weight(factor, ws[0], other.ws[0]);
sascha@1937:                 if (!Double.isNaN(wt)) {
sascha@1937:                     double q = table.getQIndex(0, km);
sascha@1937:                     if (!Double.isNaN(q)) {
sascha@1937:                         qs.add(q);
sascha@1937:                     }
sascha@1937:                 }
sascha@1937:             }
sascha@1937: 
sascha@1937:             for (int i = 1; i < ws.length; ++i) {
sascha@1937:                 double w2 = Linear.weight(factor, ws[i], other.ws[i]);
sascha@1937:                 if (Double.isNaN(w2)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937:                 if (Math.abs(w2-w) < 0.000001) {
sascha@1937:                     double q = table.getQIndex(i, km);
sascha@1937:                     if (!Double.isNaN(q)) {
sascha@1937:                         qs.add(q);
sascha@1937:                     }
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937:                 double w1 = Linear.weight(factor, ws[i-1], other.ws[i-1]);
sascha@1937:                 if (Double.isNaN(w1)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937: 
sascha@1937:                 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937: 
sascha@1937:                 double q1 = table.getQIndex(i-1, km);
sascha@1937:                 double q2 = table.getQIndex(i,   km);
sascha@1937:                 if (Double.isNaN(q1) || Double.isNaN(q2)) {
sascha@1937:                     continue;
sascha@1937:                 }
sascha@1937: 
sascha@1937:                 double q = Linear.linear(w, w1, w2, q1, q2);
sascha@1937:                 qs.add(q);
sascha@1937:             }
sascha@1937: 
sascha@1937:             return qs.toNativeArray();
sascha@1937:         }
ingo@2423: 
ingo@2423:         public double [] getMinMaxW(double [] result) {
ingo@2423:             double minW =  Double.MAX_VALUE;
ingo@2423:             double maxW = -Double.MAX_VALUE;
ingo@2423:             for (int i = 0; i < ws.length; ++i) {
ingo@2423:                 double w = ws[i];
ingo@2423:                 if (w < minW) minW = w;
ingo@2423:                 if (w > maxW) maxW = w;
ingo@2423:             }
ingo@2423:             result[0] = minW;
ingo@2423:             result[1] = maxW;
ingo@2423:             return result;
ingo@2423:         }
ingo@2423: 
ingo@2423:         public double [] getMinMaxW(Row other, double km, double [] result) {
ingo@2423:             double [] m1 = this .getMinMaxW(new double [2]);
ingo@2423:             double [] m2 = other.getMinMaxW(new double [2]);
ingo@2423:             double factor = Linear.factor(km, this.km, other.km);
ingo@2423:             result[0] = Linear.weight(factor, m1[0], m2[0]);
ingo@2423:             result[1] = Linear.weight(factor, m1[1], m2[1]);
ingo@2423:             return result;
ingo@2423:         }
sascha@633:     } // class Row
sascha@326: 
felix@1838:     /** Rows in table. */
sascha@326:     protected List<Row> rows;
sascha@326: 
felix@1838:     /** Columns in table. */
sascha@326:     protected Column [] columns;
sascha@326: 
sascha@326:     public WstValueTable() {
sascha@326:         rows = new ArrayList<Row>();
sascha@326:     }
sascha@326: 
sascha@326:     public WstValueTable(Column [] columns) {
sascha@326:         this();
sascha@326:         this.columns = columns;
sascha@326:     }
sascha@326: 
sascha@633:     public WstValueTable(Column [] columns, List<Row> rows) {
sascha@633:         this.columns = columns;
sascha@633:         this.rows    = rows;
sascha@633:     }
sascha@633: 
felix@1838:     /** Sort rows (by km). */
sascha@458:     public void sortRows() {
sascha@458:         Collections.sort(rows);
sascha@458:     }
sascha@458: 
felix@1887:     /**
felix@1887:      * @param km Given kilometer.
felix@1887:      * @param qs Given Q values.
felix@1887:      * @param ws output parameter.
felix@1887:      */
sascha@380:     public double [] interpolateW(double km, double [] qs, double [] ws) {
ingo@686:         return interpolateW(km, qs, ws, null);
ingo@686:     }
sascha@335: 
felix@1098: 
felix@1098:     /**
felix@1098:      * @param ws (output parameter), gets returned.
felix@1098:      * @return output parameter ws.
felix@1098:      */
ingo@686:     public double [] interpolateW(
ingo@686:         double      km,
sascha@742:         double []   qs,
ingo@686:         double []   ws,
ingo@686:         Calculation errors
ingo@686:     ) {
sascha@335:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@335: 
sascha@633:         QPosition qPosition = new QPosition();
sascha@633: 
sascha@335:         if (rowIndex >= 0) { // direct row match
sascha@335:             Row row = rows.get(rowIndex);
sascha@335:             for (int i = 0; i < qs.length; ++i) {
ingo@686:                 if (getQPosition(km, qs[i], qPosition) == null) {
ingo@686:                     if (errors != null) {
sascha@2166:                         errors.addProblem(km, "cannot.find.q", qs[i]);
ingo@686:                     }
ingo@686:                     ws[i] = Double.NaN;
ingo@686:                 }
ingo@686:                 else {
ingo@686:                     if (Double.isNaN(ws[i] = row.getW(qPosition))
ingo@686:                     && errors != null) {
ingo@686:                         errors.addProblem(
sascha@2166:                             km, "cannot.find.w.for.q", qs[i]);
ingo@686:                     }
ingo@686:                 }
sascha@335:             }
sascha@335:         }
sascha@335:         else { // needs bilinear interpolation
sascha@335:             rowIndex = -rowIndex -1;
sascha@335: 
sascha@335:             if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@335:                 // do not extrapolate
sascha@335:                 Arrays.fill(ws, Double.NaN);
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "km.not.found");
ingo@686:                 }
sascha@335:             }
sascha@335:             else {
sascha@633:                 Row r1 = rows.get(rowIndex-1);
sascha@633:                 Row r2 = rows.get(rowIndex);
ingo@686:                 r1.interpolateW(r2, km, qs, ws, this, errors);
sascha@335:             }
sascha@335:         }
sascha@335: 
sascha@335:         return ws;
sascha@335:     }
sascha@335: 
sascha@2559:     public double [] getMinMaxQ(double km) {
sascha@2559:         return getMinMaxQ(km, new double [2]);
sascha@2559:     }
sascha@2559: 
sascha@2559:     public double [] getMinMaxQ(double km, double [] result) {
sascha@2559:         double minQ =  Double.MAX_VALUE;
sascha@2559:         double maxQ = -Double.MAX_VALUE;
sascha@2559: 
sascha@2559:         for (int i = 0; i < columns.length; ++i) {
sascha@2559:             double q = columns[i].getQRangeTree().findQ(km);
sascha@2559:             if (!Double.isNaN(q)) {
sascha@2559:                 if (q < minQ) minQ = q;
sascha@2559:                 if (q > maxQ) maxQ = q;
sascha@2559:             }
sascha@2559:         }
sascha@2559: 
sascha@2559:         if (minQ < Double.MAX_VALUE) {
sascha@2559:             result[0] = minQ;
sascha@2559:             result[1] = maxQ;
sascha@2559:             return result;
sascha@2559:         }
sascha@2559: 
sascha@2559:         return null;
sascha@2559:     }
sascha@2559: 
sascha@2560:     public double [] getMinMaxQ(double from, double to, double step) {
sascha@2560:         double [] result = new double[2];
sascha@2560: 
sascha@2560:         double minQ =  Double.MAX_VALUE;
sascha@2560:         double maxQ = -Double.MAX_VALUE;
sascha@2560: 
sascha@2560:         if (from > to) {
sascha@2560:             double tmp = from;
sascha@2560:             from = to;
sascha@2560:             to = tmp;
sascha@2560:         }
sascha@2560: 
sascha@2560:         step = Math.max(Math.abs(step), 0.0001);
sascha@2560: 
sascha@2560:         double d = from;
sascha@2560:         for (; d <= to; d += step) {
sascha@2560:             if (getMinMaxQ(d, result) != null) {
sascha@2560:                 if (result[0] < minQ) minQ = result[0];
sascha@2560:                 if (result[1] > maxQ) maxQ = result[1];
sascha@2560:             }
sascha@2560:         }
sascha@2560: 
sascha@2560:         if (d != to) {
sascha@2560:             if (getMinMaxQ(to, result) != null) {
sascha@2560:                 if (result[0] < minQ) minQ = result[0];
sascha@2560:                 if (result[1] > maxQ) maxQ = result[1];
sascha@2560:             }
sascha@2560:         }
sascha@2560: 
sascha@2560:         return minQ < Double.MAX_VALUE
sascha@2560:             ? new double [] { minQ, maxQ }
sascha@2560:             : null;
sascha@2560:     }
sascha@2559: 
ingo@2423:     public double [] getMinMaxW(double km) {
ingo@2423:         return getMinMaxW(km, new double [2]);
ingo@2423: 
ingo@2423:     }
ingo@2423:     public double [] getMinMaxW(double km, double [] result) {
ingo@2423:         int rowIndex = Collections.binarySearch(rows, new Row(km));
ingo@2423: 
ingo@2423:         if (rowIndex >= 0) {
ingo@2423:             return rows.get(rowIndex).getMinMaxW(result);
ingo@2423:         }
ingo@2423: 
ingo@2423:         rowIndex = -rowIndex -1;
ingo@2423: 
ingo@2423:         if (rowIndex < 1 || rowIndex >= rows.size()) {
ingo@2423:             // do not extrapolate
ingo@2423:             return null;
ingo@2423:         }
ingo@2423: 
ingo@2423:         Row r1 = rows.get(rowIndex-1);
ingo@2423:         Row r2 = rows.get(rowIndex);
ingo@2423: 
ingo@2423:         return r1.getMinMaxW(r2, km, result);
ingo@2423:     }
ingo@2423: 
ingo@2423:     public double [] getMinMaxW(double from, double to, double step) {
ingo@2423:         double [] result = new double[2];
ingo@2423: 
ingo@2423:         double minW =  Double.MAX_VALUE;
ingo@2423:         double maxW = -Double.MAX_VALUE;
ingo@2423: 
ingo@2423:         if (from > to) {
ingo@2423:             double tmp = from;
ingo@2423:             from = to;
ingo@2423:             to = tmp;
ingo@2423:         }
ingo@2423: 
ingo@2423:         step = Math.max(Math.abs(step), 0.0001);
ingo@2423: 
ingo@2423:         double d = from;
ingo@2423:         for (; d <= to; d += step) {
ingo@2423:             if (getMinMaxW(d, result) != null) {
ingo@2423:                 if (result[0] < minW) minW = result[0];
ingo@2423:                 if (result[1] > maxW) maxW = result[1];
ingo@2423:             }
ingo@2423:         }
ingo@2423: 
ingo@2423:         if (d != to) {
ingo@2423:             if (getMinMaxW(to, result) != null) {
ingo@2423:                 if (result[0] < minW) minW = result[0];
ingo@2423:                 if (result[1] > maxW) maxW = result[1];
ingo@2423:             }
ingo@2423:         }
ingo@2423: 
ingo@2423:         return minW < Double.MAX_VALUE
ingo@2423:             ? new double [] { minW, maxW }
ingo@2423:             : null;
ingo@2423:     }
ingo@2423: 
felix@1838:     /**
felix@1838:      * Interpolate W and Q values at a given km.
felix@1838:      */
sascha@339:     public double [][] interpolateWQ(double km) {
ingo@686:         return interpolateWQ(km, null);
sascha@339:     }
sascha@339: 
felix@1838:     /**
felix@1838:      * Interpolate W and Q values at a given km.
felix@3265:      *
felix@3265:      * @param errors where to store errors.
felix@3265:      *
felix@3265:      * @return double double array, first index Ws, second Qs.
felix@1838:      */
ingo@686:     public double [][] interpolateWQ(double km, Calculation errors) {
ingo@686:         return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
ingo@686:     }
ingo@686: 
ingo@686: 
felix@1860:     /**
felix@1860:      * Interpolate W and Q values at a given km.
felix@1860:      */
ingo@686:     public double [][] interpolateWQ(double km, int steps, Calculation errors) {
sascha@395: 
sascha@339:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@339: 
sascha@339:         if (rowIndex >= 0) { // direct row match
sascha@339:             Row row = rows.get(rowIndex);
ingo@686:             return row.interpolateWQ(steps, this, errors);
sascha@339:         }
sascha@339: 
sascha@339:         rowIndex = -rowIndex -1;
sascha@339: 
sascha@339:         if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@339:             // do not extrapolate
ingo@686:             if (errors != null) {
sascha@2166:                 errors.addProblem(km, "km.not.found");
ingo@686:             }
sascha@339:             return new double[2][0];
sascha@339:         }
sascha@339: 
sascha@339:         Row r1 = rows.get(rowIndex-1);
sascha@339:         Row r2 = rows.get(rowIndex);
sascha@339: 
ingo@686:         return r1.interpolateWQ(r2, km, steps, this, errors);
sascha@339:     }
sascha@339: 
sascha@655:     public boolean interpolate(
sascha@655:         double    km,
sascha@655:         double [] out,
sascha@655:         QPosition qPosition,
sascha@655:         Function  qFunction
sascha@649:     ) {
sascha@427:         int R1 = rows.size()-1;
sascha@427: 
sascha@655:         out[1] = qFunction.value(getQ(qPosition, km));
sascha@655: 
sascha@655:         if (Double.isNaN(out[1])) {
sascha@655:             return false;
sascha@655:         }
sascha@633: 
sascha@633:         QPosition nPosition = new QPosition();
sascha@655:         if (getQPosition(km, out[1], nPosition) == null) {
sascha@655:             return false;
sascha@655:         }
sascha@427: 
sascha@655:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@427: 
sascha@655:         if (rowIndex >= 0) {
sascha@655:             // direct row match
sascha@655:             out[0] = rows.get(rowIndex).getW(nPosition);
sascha@655:             return !Double.isNaN(out[0]);
sascha@427:         }
sascha@655: 
sascha@655:         rowIndex = -rowIndex -1;
sascha@655: 
sascha@1681:         if (rowIndex < 1 || rowIndex > R1) {
sascha@655:             // do not extrapolate
sascha@655:             return false;
sascha@655:         }
sascha@655: 
sascha@655:         Row r1 = rows.get(rowIndex-1);
sascha@655:         Row r2 = rows.get(rowIndex);
sascha@655:         out[0] = r1.getW(r2, km, nPosition);
sascha@655: 
sascha@655:         return !Double.isNaN(out[0]);
sascha@427:     }
sascha@427: 
felix@1098: 
felix@1098:     /**
felix@1098:      * Look up interpolation of a Q at given positions.
felix@1098:      *
felix@1098:      * @param q           the non-interpolated Q value.
felix@1098:      * @param referenceKm the reference km (e.g. gauge position).
felix@1098:      * @param kms         positions for which to interpolate.
felix@1098:      * @param ws          (output) resulting interpolated ws.
felix@1098:      * @param qs          (output) resulting interpolated qs.
felix@1098:      * @param errors      calculation object to store errors.
felix@1098:      */
sascha@427:     public QPosition interpolate(
ingo@686:         double      q,
ingo@686:         double      referenceKm,
ingo@686:         double []   kms,
ingo@686:         double []   ws,
ingo@686:         double []   qs,
ingo@686:         Calculation errors
sascha@426:     ) {
ingo@686:         return interpolate(
ingo@686:             q, referenceKm, kms, ws, qs, 0, kms.length, errors);
sascha@649:     }
sascha@649: 
sascha@649:     public QPosition interpolate(
ingo@686:         double      q,
ingo@686:         double      referenceKm,
ingo@686:         double []   kms,
ingo@686:         double []   ws,
ingo@686:         double []   qs,
ingo@686:         int         startIndex,
ingo@686:         int         length,
ingo@686:         Calculation errors
sascha@649:     ) {
sascha@645:         QPosition qPosition = getQPosition(referenceKm, q);
sascha@426: 
sascha@426:         if (qPosition == null) {
sascha@633:             // we cannot locate q at km
ingo@686:             Arrays.fill(ws, Double.NaN);
ingo@686:             Arrays.fill(qs, Double.NaN);
ingo@686:             if (errors != null) {
sascha@2166:                 errors.addProblem(referenceKm, "cannot.find.q", q);
ingo@686:             }
sascha@427:             return null;
sascha@426:         }
sascha@426: 
sascha@645:         Row kmKey = new Row();
sascha@645: 
sascha@645:         int R1 = rows.size()-1;
sascha@645: 
sascha@649:         for (int i = startIndex, end = startIndex+length; i < end; ++i) {
sascha@458: 
ingo@686:             if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(kms[i], "cannot.find.q", q);
ingo@686:                 }
ingo@686:                 ws[i] = Double.NaN;
ingo@686:                 continue;
ingo@686:             }
ingo@686: 
ingo@686:             kmKey.km = kms[i];
sascha@645:             int rowIndex = Collections.binarySearch(rows, kmKey);
sascha@426: 
sascha@426:             if (rowIndex >= 0) {
sascha@426:                 // direct row match
ingo@686:                 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
ingo@686:                 && errors != null) {
sascha@2166:                     errors.addProblem(kms[i], "cannot.find.w.for.q", q);
ingo@686:                 }
sascha@426:                 continue;
sascha@426:             }
sascha@426: 
sascha@426:             rowIndex = -rowIndex -1;
sascha@426: 
sascha@1681:             if (rowIndex < 1 || rowIndex > R1) {
sascha@426:                 // do not extrapolate
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(kms[i], "km.not.found");
ingo@686:                 }
felix@1888:                 ws[i] = Double.NaN;
sascha@426:                 continue;
sascha@426:             }
sascha@426:             Row r1 = rows.get(rowIndex-1);
sascha@426:             Row r2 = rows.get(rowIndex);
sascha@633: 
ingo@686:             if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
ingo@686:             && errors != null) {
sascha@2166:                 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
ingo@686:             }
sascha@426:         }
sascha@426: 
sascha@427:         return qPosition;
sascha@426:     }
sascha@633: 
felix@1888:     /**
felix@1888:      * Linearly interpolate w at a km at a column of two rows.
felix@1888:      *
felix@1888:      * @param km   position for which to interpolate.
felix@1888:      * @param row1 first row.
felix@1888:      * @param row2 second row.
felix@1888:      * @param col  column-index at which to look.
felix@1888:      *
felix@1888:      * @return Linearly interpolated w, NaN if one of the given rows was null.
felix@1888:      */
felix@1888:     public static double linearW(double km, Row row1, Row row2, int col) {
felix@1888:         if (row1 == null || row2 == null) {
felix@1888:             return Double.NaN;
felix@1888:         }
felix@1888: 
felix@1888:         return Linear.linear(km,
felix@1888:             row1.km, row2.km,
felix@1888:             row1.ws[col], row2.ws[col]);
felix@1888:     }
felix@1888: 
felix@1888:     /**
felix@1888:      * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
felix@1888:      * of other columns).
felix@1888:      * @param km position (km) at which to interpolate/lookup.
felix@1888:      * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
felix@1888:      */
felix@1888:     public double [][] interpolateWQColumnwise(double km) {
felix@1919:         log.debug("WstValueTable.interpolateWQColumnwise");
felix@1888:         double [] qs = new double[columns.length];
felix@1888:         double [] ws = new double[columns.length];
felix@1888: 
felix@1888:         // Find out row from where we will start searching.
felix@1888:         int rowIndex = Collections.binarySearch(rows, new Row(km));
felix@1888: 
felix@1888:         if (rowIndex < 0) {
felix@1888:             rowIndex = -rowIndex -1;
felix@1888:         }
felix@2248: 
felix@2248:         // TODO Beyond definition, we could stop more clever.
sascha@1937:         if (rowIndex >= rows.size()) {
felix@1919:             rowIndex = rows.size() -1;
sascha@1937:         }
felix@1888: 
felix@1888:         Row startRow = rows.get(rowIndex);
felix@1888: 
felix@1888:         for (int col = 0; col < columns.length; col++) {
felix@1888:             qs[col] = columns[col].getQRangeTree().findQ(km);
felix@1888:             if (startRow.km == km && startRow.ws[col] != Double.NaN) {
felix@1888:                 // Great. W is defined at km.
felix@1888:                 ws[col] = startRow.ws[col];
felix@1888:                 continue;
felix@1888:             }
felix@1888: 
felix@1888:             // Search neighbouring rows that define w at this col.
felix@1888:             Row rowBefore = null;
felix@1888:             Row rowAfter  = null;
felix@1899:             for (int before = rowIndex -1; before >= 0; before--) {
felix@1888:                 if (!Double.isNaN(rows.get(before).ws[col])) {
felix@1888:                     rowBefore = rows.get(before);
felix@1888:                     break;
felix@1888:                 }
felix@1888:             }
felix@2248:             if (rowBefore != null) {
felix@2248:                 for (int after = rowIndex, R = rows.size(); after < R; after++) {
felix@2248:                     if (!Double.isNaN(rows.get(after).ws[col])) {
felix@2248:                         rowAfter = rows.get(after);
felix@2248:                         break;
felix@2248:                     }
felix@1888:                 }
felix@1888:             }
felix@1888: 
felix@1888:             ws[col] = linearW(km, rowBefore, rowAfter, col);
felix@1888:         }
felix@1888: 
felix@1888:         return new double [][] {qs, ws};
felix@1888:     }
felix@1888: 
sascha@1937:     public double [] findQsForW(double km, double w) {
sascha@1937: 
sascha@1937:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@1937: 
sascha@1937:         if (rowIndex >= 0) {
sascha@1937:             return rows.get(rowIndex).findQsForW(w, this);
sascha@1937:         }
sascha@1937: 
sascha@1937:         rowIndex = -rowIndex - 1;
sascha@1937: 
sascha@1937:         if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@1939:             // Do not extrapolate.
sascha@1937:             return new double[0];
sascha@1937:         }
sascha@1937: 
sascha@1939:         // Needs bilinear interpolation.
sascha@1937:         Row r1 = rows.get(rowIndex-1);
sascha@1937:         Row r2 = rows.get(rowIndex);
sascha@1937: 
sascha@1937:         return r1.findQsForW(r2, w, km, this);
sascha@1937:     }
sascha@1937: 
sascha@1939:     protected SplineFunction createSpline(double km, Calculation errors) {
sascha@1939: 
sascha@1939:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@1939: 
sascha@1939:         if (rowIndex >= 0) {
sascha@1939:             SplineFunction sf = rows.get(rowIndex).createSpline(this, errors);
sascha@1939:             if (sf == null && errors != null) {
sascha@2166:                 errors.addProblem(km, "cannot.create.wq.relation");
sascha@1939:             }
sascha@1939:             return sf;
sascha@1939:         }
sascha@1939: 
sascha@1939:         rowIndex = -rowIndex - 1;
sascha@1939: 
sascha@1939:         if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@1939:             // Do not extrapolate.
sascha@1939:             if (errors != null) {
sascha@2166:                 errors.addProblem(km, "km.not.found");
sascha@1939:             }
sascha@1939:             return null;
sascha@1939:         }
sascha@1939: 
sascha@1939:         // Needs bilinear interpolation.
sascha@1939:         Row r1 = rows.get(rowIndex-1);
sascha@1939:         Row r2 = rows.get(rowIndex);
sascha@1939: 
sascha@1939:         SplineFunction sf = r1.createSpline(r2, km, this, errors);
sascha@1939:         if (sf == null && errors != null) {
sascha@2166:             errors.addProblem(km, "cannot.create.wq.relation");
sascha@1939:         }
sascha@1939: 
sascha@1939:         return sf;
sascha@1939:     }
sascha@1939: 
sascha@1939:     /** 'Bezugslinienverfahren' */
sascha@1939:     public double [][] relateWs(
sascha@3076:         double      km1,
sascha@1939:         double      km2,
sascha@2186:         Calculation errors
sascha@2186:     ) {
sascha@2186:         return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
sascha@2186:     }
sascha@2186: 
ingo@2330:     private static class ErrorHandler {
ingo@2330: 
ingo@2330:         boolean     hasErrors;
ingo@2330:         Calculation errors;
ingo@2330: 
ingo@2330:         ErrorHandler(Calculation errors) {
ingo@2330:             this.errors = errors;
ingo@2330:         }
ingo@2330: 
ingo@2330:         void error(double km, String key, Object ... args) {
ingo@2330:             if (errors != null && !hasErrors) {
ingo@2330:                 hasErrors = true;
ingo@2330:                 errors.addProblem(km, key, args);
ingo@2330:             }
ingo@2330:         }
ingo@2330:     } // class ErrorHandler
ingo@2330: 
ingo@2330: 
sascha@2186:     /* TODO: Add optimized methods of relateWs to relate one
sascha@3076:      *       start km to many end kms. The index generation/spline stuff for
sascha@2186:      *       the start km is always the same.
sascha@2186:      */
sascha@2186:     public double [][] relateWs(
sascha@3076:         double      km1,
sascha@2186:         double      km2,
sascha@1939:         int         numSamples,
sascha@1939:         Calculation errors
sascha@1939:     ) {
sascha@1939:         SplineFunction sf1 = createSpline(km1, errors);
sascha@1939:         if (sf1 == null) {
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
felix@2248:         SplineFunction sf2 = createSpline(km2, errors);
sascha@1939:         if (sf2 == null) {
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
sascha@1939:         PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
sascha@1939:         if (iQ1 == null) {
sascha@1939:             if (errors != null) {
sascha@2166:                 errors.addProblem(km1, "cannot.create.index.q.relation");
sascha@1939:             }
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
sascha@2404:         PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
sascha@1939:         if (iQ2 == null) {
sascha@1939:             if (errors != null) {
sascha@2166:                 errors.addProblem(km2, "cannot.create.index.q.relation");
sascha@1939:             }
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
sascha@2316:         int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
sascha@1939:         double stepWidth = N/(double)numSamples;
sascha@1939: 
sascha@1939:         PolynomialSplineFunction qW1 = sf1.spline;
sascha@1939:         PolynomialSplineFunction qW2 = sf2.spline;
sascha@1939: 
sascha@2316:         TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
sascha@2316:         TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
sascha@2316:         TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
sascha@2316:         TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);
sascha@1939: 
ingo@2330:         ErrorHandler err = new ErrorHandler(errors);
sascha@1939: 
sascha@2316:         int i = 0;
sascha@2316:         for (double p = 0d; p <= N-1; p += stepWidth, ++i) {
ingo@2330: 
ingo@2330:             double q1;
sascha@1939:             try {
ingo@2330:                 q1 = iQ1.value(p);
sascha@1939:             }
sascha@1939:             catch (ArgumentOutsideDomainException aode) {
ingo@2330:                 err.error(km1, "w.w.qkm1.failed", p);
ingo@2330:                 continue;
sascha@1939:             }
ingo@2330: 
ingo@2330:             double w1;
ingo@2330:             try {
ingo@2330:                 w1 = qW1.value(q1);
ingo@2330:             }
ingo@2330:             catch (ArgumentOutsideDomainException aode) {
sascha@2403:                 err.error(km1, "w.w.wkm1.failed", q1, p);
ingo@2330:                 continue;
ingo@2330:             }
ingo@2330: 
ingo@2330:             double q2;
ingo@2330:             try {
ingo@2330:                 q2 = iQ2.value(p);
ingo@2330:             }
ingo@2330:             catch (ArgumentOutsideDomainException aode) {
ingo@2330:                 err.error(km2, "w.w.qkm2.failed", p);
ingo@2330:                 continue;
ingo@2330:             }
ingo@2330: 
ingo@2330:             double w2;
ingo@2330:             try {
ingo@2330:                 w2 = qW2.value(q2);
ingo@2330:             }
ingo@2330:             catch (ArgumentOutsideDomainException aode) {
sascha@2403:                 err.error(km2, "w.w.wkm2.failed", q2, p);
ingo@2330:                 continue;
ingo@2330:             }
ingo@2330: 
ingo@2330:             ws1.add(w1);
ingo@2330:             ws2.add(w2);
ingo@2330:             qs1.add(q1);
ingo@2330:             qs2.add(q2);
sascha@1939:         }
sascha@1939: 
sascha@2316:         return new double [][] {
sascha@2316:             ws1.toNativeArray(),
sascha@3076:             qs1.toNativeArray(),
sascha@2316:             ws2.toNativeArray(),
sascha@2316:             qs2.toNativeArray() };
sascha@1939:     }
sascha@1939: 
sascha@633:     public QPosition getQPosition(double km, double q) {
sascha@633:         return getQPosition(km, q, new QPosition());
sascha@633:     }
sascha@633: 
sascha@633:     public QPosition getQPosition(double km, double q, QPosition qPosition) {
sascha@633: 
sascha@633:         if (columns.length == 0) {
sascha@633:             return null;
sascha@633:         }
sascha@633: 
sascha@633:         double qLast = columns[0].getQRangeTree().findQ(km);
sascha@633: 
sascha@633:         if (Math.abs(qLast - q) < 0.00001) {
sascha@633:             return qPosition.set(0, 1d);
sascha@633:         }
sascha@633: 
sascha@633:         for (int i = 1; i < columns.length; ++i) {
sascha@633:             double qCurrent = columns[i].getQRangeTree().findQ(km);
sascha@633:             if (Math.abs(qCurrent - q) < 0.00001) {
sascha@633:                 return qPosition.set(i, 1d);
sascha@633:             }
sascha@633: 
sascha@633:             double qMin, qMax;
sascha@633:             if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
sascha@633:             else                  { qMin = qCurrent; qMax = qLast; }
sascha@633: 
sascha@633:             if (q > qMin && q < qMax) {
sascha@655:                 double weight = Linear.factor(q, qLast, qCurrent);
sascha@633:                 return qPosition.set(i, weight);
sascha@633:             }
sascha@633:             qLast = qCurrent;
sascha@633:         }
sascha@633: 
sascha@633:         return null;
sascha@633:     }
sascha@633: 
sascha@633:     public double getQIndex(int index, double km) {
sascha@633:         return columns[index].getQRangeTree().findQ(km);
sascha@633:     }
sascha@633: 
sascha@633:     public double getQ(QPosition qPosition, double km) {
sascha@633:         int    index  = qPosition.index;
sascha@633:         double weight = qPosition.weight;
sascha@633: 
sascha@633:         if (weight == 1d) {
sascha@633:             return columns[index].getQRangeTree().findQ(km);
sascha@633:         }
sascha@633:         double q1 = columns[index-1].getQRangeTree().findQ(km);
sascha@633:         double q2 = columns[index  ].getQRangeTree().findQ(km);
sascha@655:         return Linear.weight(weight, q1, q2);
sascha@633:     }
sascha@3736: 
sascha@3736:     public double [][] interpolateTabulated(double km) {
sascha@3736:         return interpolateTabulated(km, new double[2][columns.length]);
sascha@3736:     }
sascha@3736: 
sascha@3736:     public double [][] interpolateTabulated(double km, double [][] result) {
sascha@3736: 
sascha@3736:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@3736: 
sascha@3736:         if (rowIndex >= 0) {
sascha@3736:             // Direct hit -> copy ws.
sascha@3736:             Row row = rows.get(rowIndex);
sascha@3736:             System.arraycopy(
sascha@3736:                 row.ws, 0, result[0], 0,
sascha@3736:                 Math.min(row.ws.length, result[0].length));
sascha@3736:         }
sascha@3736:         else {
sascha@3736:             rowIndex = -rowIndex -1;
sascha@3736:             if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@3736:                 // Out of bounds.
sascha@3736:                 return null;
sascha@3736:             }
sascha@3736:             // Interpolate ws.
sascha@3736:             Row r1 = rows.get(rowIndex-1);
sascha@3736:             Row r2 = rows.get(rowIndex);
sascha@3736:             double factor = Linear.factor(km, r1.km, r2.km);
sascha@3736:             Linear.weight(factor, r1.ws, r2.ws, result[0]);
sascha@3736:         }
sascha@3736: 
sascha@3736:         double [] qs = result[1];
sascha@3736:         for (int i = Math.min(qs.length, columns.length)-1; i >= 0; --i) {
sascha@3736:             qs[i] = columns[i].getQRangeTree().findQ(km);
sascha@3736:         }
sascha@3736:         return result;
sascha@3736:     }
sascha@3743: 
felix@4119: 
felix@4119:     /** Find ranges that are between km1 and km2 (inclusive?) */
sascha@3743:     public List<Range> findSegments(double km1, double km2) {
sascha@3743:         return columns.length != 0
sascha@3743:             ? columns[columns.length-1].getQRangeTree().findSegments(km1, km2)
sascha@3743:             : Collections.<Range>emptyList();
sascha@3743:     }
sascha@326: }
sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :