Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java @ 3318:dbe2f85bf160
merged flys-artifacts/2.8
author | Thomas Arendsen Hein <thomas@intevation.de> |
---|---|
date | Fri, 28 Sep 2012 12:14:35 +0200 |
parents | a0d9a99a5d17 |
children | cb11919cccf9 |
line wrap: on
line source
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 :