--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java	Fri Sep 28 12:14:35 2012 +0200
@@ -0,0 +1,322 @@
+package de.intevation.flys.artifacts.model;
+
+import java.util.List;
+import java.util.Arrays;
+import java.util.Comparator;
+import java.util.Collections;
+
+import de.intevation.flys.utils.DoubleUtil;
+
+import de.intevation.flys.model.River;
+import de.intevation.flys.model.Gauge;
+import de.intevation.flys.model.DischargeTable;
+
+import de.intevation.flys.artifacts.model.WstValueTable.QPosition;
+
+import de.intevation.flys.artifacts.math.Function;
+import de.intevation.flys.artifacts.math.Linear;
+import de.intevation.flys.artifacts.math.Identity;
+import de.intevation.flys.artifacts.math.BackJumpCorrector;
+
+import org.apache.log4j.Logger;
+
+public class Calculation4
+extends      Calculation
+{
+    private static Logger logger = Logger.getLogger(Calculation4.class);
+
+    public static final double MINIMAL_STEP_WIDTH = 1e-5;
+
+    public static final Comparator<Segment> REF_CMP =
+        new Comparator<Segment>() {
+            public int compare(Segment a, Segment b) {
+                double d = a.referencePoint - b.referencePoint;
+                if (d < 0d) return -1;
+                return d > 0d ? +1 : 0;
+            }
+        };
+
+    protected List<Segment> segments;
+
+    protected boolean       isQ;
+
+    public Calculation4() {
+    }
+
+    public Calculation4(List<Segment> segments, River river, boolean isQ) {
+
+        this.segments = segments;
+        this.isQ      = isQ;
+
+        int numResults = -1;
+
+        // assign reference points
+        for (Segment segment: segments) {
+            Gauge gauge = river.maxOverlap(segment.getFrom(), segment.getTo());
+
+            if (gauge == null) {
+                logger.warn("no gauge found. Defaults to mid point.");
+                segment.setReferencePoint(
+                    0.5*(segment.getFrom()+segment.getTo()));
+            }
+            else {
+                double ref = gauge.getStation().doubleValue();
+                logger.debug(
+                    "reference gauge: " + gauge.getName() +
+                    " (km " + ref + ")");
+                segment.setReferencePoint(ref);
+            }
+
+            double [] values = segment.values;
+
+            if (numResults == -1) {
+                numResults = values.length;
+            }
+            else if (numResults != values.length) {
+                throw new IllegalArgumentException("wrong length of values");
+            }
+
+            // convert to Q if needed
+            if (!isQ && gauge != null) {
+
+                DischargeTable dt = gauge.fetchMasterDischargeTable();
+
+                double [][] table =
+                    DischargeTables.loadDischargeTableValues(dt, 1);
+
+                // need the original values for naming
+                segment.backup();
+
+                for (int i = 0; i < values.length; ++i) {
+                    double w = values[i] / 100.0;
+                    double [] qs = DischargeTables.getQsForW(table, w);
+                    if (qs.length == 0) {
+                        logger.warn("No Qs found for W = " + values[i]);
+                        addProblem("cannot.find.w.for.q", values[i]);
+                        values[i] = Double.NaN;
+                    }
+                    else {
+                        values[i] = qs[0];
+                        if (qs.length > 1) {
+                            logger.warn(
+                                "More than one Q found for W = " + values[i]);
+                        }
+                    }
+                }
+            }
+        } // for all segments
+
+        Collections.sort(segments, REF_CMP);
+    }
+
+    public CalculationResult calculate(
+        WstValueTable table,
+        double from, double to, double step
+    ) {
+        boolean debug = logger.isDebugEnabled();
+
+        if (debug) {
+            logger.debug(
+                "calculate from " + from + " to " + to + " step " + step);
+            logger.debug("# segments: " + segments.size());
+            for (Segment segment: segments) {
+                logger.debug("  " + segment);
+            }
+        }
+
+        if (segments.isEmpty()) {
+            logger.debug("no segments found");
+            addProblem("no.segments.found");
+            return new CalculationResult(new WQKms[0], this);
+        }
+
+        int numResults = segments.get(0).values.length;
+
+        if (numResults < 1) {
+            logger.debug("no values given");
+            addProblem("no.values.given");
+            return new CalculationResult(new WQKms[0], this);
+        }
+
+
+        WQKms [] results = new WQKms[numResults];
+        for (int i = 0; i < results.length; ++i) {
+            results[i] = new WQKms();
+        }
+
+        if (Math.abs(step) < MINIMAL_STEP_WIDTH) {
+            step = MINIMAL_STEP_WIDTH;
+        }
+
+        if (from > to) {
+            step = -step;
+        }
+
+        QPosition [] qPositions = new QPosition[numResults];
+
+        Function [] functions = new Function[numResults];
+
+        double [] out = new double[2];
+
+        Segment sentinel = new Segment(Double.MAX_VALUE);
+        Segment s1 = sentinel, s2 = sentinel;
+
+        for (double pos = from;
+             from < to ? pos <= to : pos >= to;
+             pos = DoubleUtil.round(pos + step)
+        ) {
+            if (pos < s1.referencePoint || pos > s2.referencePoint) {
+                if (debug) {
+                    logger.debug("need to find new interval for " + pos);
+                }
+                // find new interval
+                if (pos <= segments.get(0).referencePoint) {
+                    // before first segment -> "gleichwertig"
+                    if (debug) {
+                        logger.debug("before first segment -> gleichwertig");
+                    }
+                    Segment   first  = segments.get(0);
+                    double [] values = first.values;
+                    double    refPos = first.referencePoint;
+                    for (int i = 0; i < qPositions.length; ++i) {
+                        qPositions[i] = table.getQPosition(
+                            refPos, values[i]);
+                    }
+                    sentinel.setReferencePoint(-Double.MAX_VALUE);
+                    s1 = sentinel;
+                    s2 = segments.get(0);
+                    Arrays.fill(functions, Identity.IDENTITY);
+                }
+                else if (pos >= segments.get(segments.size()-1).referencePoint) {
+                    // after last segment -> "gleichwertig"
+                    if (debug) {
+                        logger.debug("after last segment -> gleichwertig");
+                    }
+                    Segment   last   = segments.get(segments.size()-1);
+                    double [] values = last.values;
+                    double    refPos = last.referencePoint;
+                    for (int i = 0; i < qPositions.length; ++i) {
+                        qPositions[i] = table.getQPosition(
+                            refPos, values[i]);
+                    }
+                    sentinel.setReferencePoint(Double.MAX_VALUE);
+                    s1 = last;
+                    s2 = sentinel;
+                    Arrays.fill(functions, Identity.IDENTITY);
+                }
+                else { // "ungleichwertig"
+                    // find matching interval
+                    if (debug) {
+                        logger.debug("in segments -> ungleichwertig");
+                    }
+                    s1 = s2 = null;
+                    for (int i = 1, N = segments.size(); i < N; ++i) {
+                        Segment si1 = segments.get(i-1);
+                        Segment si  = segments.get(i);
+                        if (debug) {
+                            logger.debug("check " + pos + " in " +
+                                si1.referencePoint + " - " + si.referencePoint);
+                        }
+                        if (pos >= si1.referencePoint
+                        &&  pos <= si. referencePoint) {
+                            s1 = si1;
+                            s2 = si;
+                            break;
+                        }
+                    }
+
+                    if (s1 == null) {
+                        throw new IllegalStateException("no interval found");
+                    }
+
+                    Segment anchor, free;
+
+                    if (from > to) { anchor = s1; free = s2; }
+                    else           { anchor = s2; free = s1; }
+
+                    // build transforms based on "gleichwertiger" phase
+                    for (int i = 0; i < qPositions.length; ++i) {
+                        QPosition qi = table.getQPosition(
+                            anchor.referencePoint,
+                            anchor.values[i]);
+
+                        if ((qPositions[i] = qi) == null) {
+                            addProblem(pos, "cannot.find.q", anchor.values[i]);
+                            functions[i] = Identity.IDENTITY;
+                        }
+                        else {
+                            double qA = table.getQ(qi, anchor.referencePoint);
+                            double qF = table.getQ(qi, free  .referencePoint);
+
+                            functions[i] = Double.isNaN(qA) || Double.isNaN(qF)
+                                ? Identity.IDENTITY
+                                : new Linear(
+                                    qA, qF,
+                                    anchor.values[i], free.values[i]);
+
+                            if (debug) {
+                                logger.debug(
+                                    anchor.referencePoint + ": " +
+                                    qA + " -> " + functions[i].value(qA) +
+                                    " / " + free.referencePoint + ": " +
+                                    qF + " -> " + functions[i].value(qF));
+                            }
+                        }
+                    } // build transforms
+                } // "ungleichwertiges" interval
+            } // find matching interval
+
+            for (int i = 0; i < qPositions.length; ++i) {
+                QPosition qPosition = qPositions[i];
+
+                if (qPosition == null) {
+                    continue;
+                }
+
+                if (table.interpolate(pos, out, qPosition, functions[i])) {
+                    results[i].add(out[0], out[1], pos);
+                }
+                else {
+                    addProblem(pos, "cannot.interpolate.w.q");
+                }
+            }
+        }
+
+        // Backjump correction
+        for (int i = 0; i < results.length; ++i) {
+            BackJumpCorrector bjc = new BackJumpCorrector();
+
+            double [] ws  = results[i].getWs();
+            double [] kms = results[i].getKms();
+
+            if (bjc.doCorrection(kms, ws, this)) {
+                results[i] = new WQCKms(results[i], bjc.getCorrected());
+            }
+        }
+
+        // name the curves
+        for (int i = 0; i < results.length; ++i) {
+            results[i].setName(createName(i));
+        }
+
+        return new CalculationResult(results, this);
+    }
+
+    protected String createName(int index) {
+        // TODO: I18N
+        StringBuilder sb = new StringBuilder(isQ ? "Q" : "W");
+        sb.append(" benutzerdefiniert (");
+        for (int i = 0, N = segments.size(); i < N; ++i) {
+            if (i > 0) {
+                sb.append("; ");
+            }
+            Segment segment = segments.get(i);
+            sb.append((segment.backup != null
+                ? segment.backup
+                : segment.values)[index]);
+        }
+        sb.append(')');
+        return sb.toString();
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :