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