teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde teichmann@5863: * Software engineering by Intevation GmbH teichmann@5863: * teichmann@5994: * This file is Free Software under the GNU AGPL (>=v3) teichmann@5863: * and comes with ABSOLUTELY NO WARRANTY! Check out the teichmann@5994: * documentation coming with Dive4Elements River for details. teichmann@5863: */ teichmann@5863: teichmann@5831: package org.dive4elements.river.exports.fixings; sascha@3220: teichmann@5831: import org.dive4elements.artifacts.CallMeta; sascha@3220: teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.Function; sascha@3220: teichmann@5831: import org.dive4elements.river.artifacts.model.Parameters; sascha@3220: teichmann@5831: import org.dive4elements.river.artifacts.resources.Resources; teichmann@5831: teichmann@5831: import org.dive4elements.river.exports.ATWriter; sascha@3220: aheinecke@6325: import org.dive4elements.river.model.Gauge; aheinecke@6325: import org.dive4elements.river.model.River; aheinecke@6325: sascha@3220: import java.io.IOException; sascha@3220: import java.io.PrintWriter; sascha@3220: import java.io.Writer; sascha@3220: sascha@3220: import java.util.Locale; aheinecke@7687: import java.util.Arrays; sascha@3220: sascha@3220: import org.apache.log4j.Logger; sascha@3220: felix@3905: /** Export Fixation Analysis Results to AT. */ sascha@3220: public class FixATWriter sascha@3220: { teichmann@8202: /** Private log. */ sascha@3220: private static Logger log = Logger.getLogger(FixATWriter.class); sascha@3220: sascha@3220: public static final String I18N_HEADER_KEY = sascha@3220: "fix.export.at.header"; sascha@3220: aheinecke@6325: public static final String I18N_GAUGE_HEADER_KEY = aheinecke@6325: "fix.export.at.gauge.header"; aheinecke@6325: sascha@3220: public static final String I18N_HEADER_DEFAULT = sascha@3220: "Exported fixings discharge curve for {0} {0}-km: {1}"; sascha@3220: aheinecke@6325: public static final String I18N_GAUGE_HEADER_DEFAULT = aheinecke@6325: "Exported fixings discharge curve for {0}, gauge: {1} datum[{3}] = {2}"; aheinecke@6325: sascha@3220: public static final String [] Q_MAX_COLUMN = new String [] { "max_q" }; sascha@3220: sascha@3641: private static final int MAX_ITERATIONS = 10000; sascha@3641: private static final double EPSILON = 1e-8; sascha@3643: private static final double MIN_Q = 1e-4; sascha@3220: sascha@3220: protected Function function; sascha@3220: protected Parameters parameters; sascha@3220: sascha@3220: public FixATWriter() { sascha@3220: } sascha@3220: sascha@3220: public FixATWriter(Function function, Parameters parameters) { sascha@3220: this.function = function; sascha@3220: this.parameters = parameters; sascha@3220: } sascha@3220: sascha@3220: public void write( sascha@3220: Writer writer, sascha@3220: CallMeta meta, aheinecke@6325: River river, sascha@3220: double km sascha@3220: ) sascha@3220: throws IOException { sascha@3220: PrintWriter out = new PrintWriter(writer); aheinecke@6325: aheinecke@6325: int subtractPNP = 0; aheinecke@6325: // Special case handling for at's at gauges aheinecke@6325: Gauge gauge = river.determineGaugeByPosition(km); aheinecke@6325: if (Math.abs(km - gauge.getStation().doubleValue()) < 1e-4) { aheinecke@6325: printGaugeHeader(out, meta, river, gauge); aheinecke@6325: subtractPNP = (int)Math.round(gauge.getDatum().doubleValue() * 100); aheinecke@6325: } else { aheinecke@6325: printHeader(out, meta, river.getName(), km); aheinecke@6325: } sascha@3220: sascha@3220: double [] coeffs = parameters.interpolate( sascha@3220: "km", km, function.getParameterNames()); sascha@3220: sascha@3220: double [] qMax = parameters.interpolate( sascha@3220: "km", km, Q_MAX_COLUMN); sascha@3220: sascha@3220: if (coeffs == null || qMax == null) { sascha@3220: log.debug("No data found at km " + km + "."); sascha@3220: return; sascha@3220: } sascha@3220: teichmann@5831: org.dive4elements.river.artifacts.math.Function funcInst = sascha@3220: function.instantiate(coeffs); sascha@3220: sascha@3642: // Increase Q max about 5%. sascha@3642: qMax[0] += Math.abs(qMax[0])*0.05; sascha@3642: sascha@3220: double wMax = funcInst.value(qMax[0]); sascha@3220: sascha@3220: if (Double.isNaN(wMax) || wMax < 0d) { sascha@3220: log.debug("function '" + function.getName() + sascha@3220: "' eval failed at " + wMax); sascha@3220: return; sascha@3220: } sascha@3220: sascha@3220: Function inverse = function.getInverse(); sascha@3220: teichmann@5831: org.dive4elements.river.artifacts.math.Function invInst = sascha@3220: inverse.instantiate(coeffs); sascha@3220: sascha@3641: double wMin = minW(invInst, wMax, qMax[0]); sascha@3220: aheinecke@6637: double wMinCM = wMin * 100d - subtractPNP; sascha@3641: double wMaxCM = wMax * 100d; sascha@3220: sascha@3641: if ((wMinCM - (int)wMinCM) > 0d) { sascha@3641: wMinCM = (int)wMinCM + 1d; sascha@3641: } sascha@3641: aheinecke@7687: int wRow = ((int)wMinCM / 10) * 10; aheinecke@7687: aheinecke@6637: double w = (wMinCM + subtractPNP) / 100.0; sascha@3641: sascha@3641: int wcm = ((int)wMinCM) % 10; sascha@3641: sascha@3641: if (log.isDebugEnabled()) { aheinecke@7687: log.debug("km: " + km); sascha@3641: log.debug("wMinCM: " + wMinCM); sascha@3641: log.debug("wMaxCM: " + wMaxCM); sascha@3641: log.debug("wcm: " + wcm); aheinecke@6637: log.debug("subtractPNP: " + subtractPNP); aheinecke@7687: log.debug("coeffs: " + Arrays.toString(coeffs)); aheinecke@7687: log.debug("function description: " + inverse.getDescription()); sascha@3641: } sascha@3641: aheinecke@6637: out.printf(Locale.US, "%8d", wRow); sascha@3641: aheinecke@6209: for (int i = 0; i < wcm; i++) { aheinecke@6209: out.print(ATWriter.EMPTY); sascha@3220: } sascha@3220: sascha@3641: for (;;) { sascha@3641: while (wcm++ < 10) { sascha@3641: if (w > wMax) { sascha@3641: break; sascha@3641: } sascha@3641: double q = invInst.value(w); sascha@3641: if (Double.isNaN(w)) { sascha@3641: out.print(ATWriter.EMPTY); sascha@3641: } sascha@3641: else { sascha@3641: ATWriter.printQ(out, q); sascha@3641: } sascha@3641: w += 0.01d; sascha@3641: } sascha@3641: out.println(); sascha@3641: if (w > wMax) { sascha@3641: break; sascha@3641: } aheinecke@6637: out.printf(Locale.US, "%8d", (wRow += 10)); sascha@3641: wcm = 0; sascha@3220: } sascha@3220: sascha@3220: out.flush(); sascha@3220: } sascha@3220: sascha@3220: protected void printHeader( sascha@3220: PrintWriter out, sascha@3220: CallMeta meta, sascha@3220: String river, sascha@3220: double km sascha@3220: ) { aheinecke@6208: out.println("*" + Resources.format( sascha@3220: meta, sascha@3220: I18N_HEADER_KEY, sascha@3220: I18N_HEADER_DEFAULT, sascha@3220: river, km)); sascha@3220: } sascha@3641: aheinecke@6325: protected void printGaugeHeader( aheinecke@6325: PrintWriter out, aheinecke@6325: CallMeta meta, aheinecke@6325: River river, aheinecke@6325: Gauge gauge aheinecke@6325: ) { aheinecke@6325: out.println("*" + Resources.format( aheinecke@6325: meta, aheinecke@6325: I18N_GAUGE_HEADER_KEY, aheinecke@6325: I18N_GAUGE_HEADER_DEFAULT, aheinecke@6325: new Object[] { river.getName(), gauge.getName(), aheinecke@6325: gauge.getDatum(), river.getWstUnit().getName() })); aheinecke@6325: } aheinecke@6325: sascha@3641: private static double minW( teichmann@5831: org.dive4elements.river.artifacts.math.Function function, sascha@3641: double maxW, sascha@3641: double maxQ sascha@3641: ) { sascha@3641: double stepWidth = 10d; sascha@3641: sascha@3641: double lastW = maxW; sascha@3641: double lastQ = maxQ; sascha@3641: sascha@3641: for (int i = 0; i < MAX_ITERATIONS; ++i) { sascha@3641: double w = lastW - stepWidth; aheinecke@7687: aheinecke@7687: if (w <= 0) { aheinecke@7687: return 0; aheinecke@7687: } aheinecke@7687: sascha@3641: double q = function.value(w); sascha@3641: sascha@3643: if (Double.isNaN(q) || q > lastQ || q < MIN_Q) { sascha@3641: if (stepWidth < EPSILON) { sascha@3641: break; sascha@3641: } sascha@3641: stepWidth *= 0.5d; sascha@3641: continue; sascha@3641: } sascha@3641: sascha@3641: lastW = w; sascha@3641: lastQ = q; sascha@3641: } sascha@3641: sascha@3641: return lastW; sascha@3641: } sascha@3220: } sascha@3220: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :