Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/exports/fixings/FixATWriter.java @ 3643:171db4d5d3cb
FixA: AT export: prevent some numerical problems with steep functions around zero.
flys-artifacts/trunk@5363 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 05 Sep 2012 08:38:12 +0000 |
parents | 83c0735092a3 |
children | 8bbb9e173297 |
line wrap: on
line source
package de.intevation.flys.exports.fixings; import de.intevation.artifacts.CallMeta; import de.intevation.flys.artifacts.math.fitting.Function; import de.intevation.flys.artifacts.model.Parameters; import de.intevation.flys.artifacts.resources.Resources; import de.intevation.flys.exports.ATWriter; import java.io.IOException; import java.io.PrintWriter; import java.io.Writer; import java.util.Locale; import org.apache.log4j.Logger; public class FixATWriter { private static Logger log = Logger.getLogger(FixATWriter.class); public static final String I18N_HEADER_KEY = "fix.export.at.header"; public static final String I18N_HEADER_DEFAULT = "Exported fixings discharge curve for {0} {0}-km: {1}"; public static final String [] Q_MAX_COLUMN = new String [] { "max_q" }; private static final int MAX_ITERATIONS = 10000; private static final double EPSILON = 1e-8; private static final double MIN_Q = 1e-4; protected Function function; protected Parameters parameters; public FixATWriter() { } public FixATWriter(Function function, Parameters parameters) { this.function = function; this.parameters = parameters; } public void write( Writer writer, CallMeta meta, String river, double km ) throws IOException { PrintWriter out = new PrintWriter(writer); printHeader(out, meta, river, km); double [] coeffs = parameters.interpolate( "km", km, function.getParameterNames()); double [] qMax = parameters.interpolate( "km", km, Q_MAX_COLUMN); if (coeffs == null || qMax == null) { log.debug("No data found at km " + km + "."); return; } de.intevation.flys.artifacts.math.Function funcInst = function.instantiate(coeffs); // Increase Q max about 5%. qMax[0] += Math.abs(qMax[0])*0.05; double wMax = funcInst.value(qMax[0]); if (Double.isNaN(wMax) || wMax < 0d) { log.debug("function '" + function.getName() + "' eval failed at " + wMax); return; } Function inverse = function.getInverse(); de.intevation.flys.artifacts.math.Function invInst = inverse.instantiate(coeffs); double wMin = minW(invInst, wMax, qMax[0]); double wMinCM = wMin * 100d; double wMaxCM = wMax * 100d; int wRow = ((int)wMinCM / 10) * 10; if ((wMinCM - (int)wMinCM) > 0d) { wMinCM = (int)wMinCM + 1d; } double w = wMinCM / 100.0; int wcm = ((int)wMinCM) % 10; if (log.isDebugEnabled()) { log.debug("wMinCM: " + wMinCM); log.debug("wMaxCM: " + wMaxCM); log.debug("wcm: " + wcm); } out.printf(Locale.US, "%8d", wRow); if (wcm > 0) { int rest = 10 - wcm; while (rest-- > 0) { out.print(ATWriter.EMPTY); } } for (;;) { while (wcm++ < 10) { if (w > wMax) { break; } double q = invInst.value(w); if (Double.isNaN(w)) { out.print(ATWriter.EMPTY); } else { ATWriter.printQ(out, q); } w += 0.01d; } out.println(); if (w > wMax) { break; } out.printf(Locale.US, "%8d", wRow += 10); wcm = 0; } out.flush(); } protected void printHeader( PrintWriter out, CallMeta meta, String river, double km ) { out.println(Resources.format( meta, I18N_HEADER_KEY, I18N_HEADER_DEFAULT, river, km)); } private static double minW( de.intevation.flys.artifacts.math.Function function, double maxW, double maxQ ) { double stepWidth = 10d; double lastW = maxW; double lastQ = maxQ; for (int i = 0; i < MAX_ITERATIONS; ++i) { double w = lastW - stepWidth; double q = function.value(w); if (Double.isNaN(q) || q > lastQ || q < MIN_Q) { if (stepWidth < EPSILON) { break; } stepWidth *= 0.5d; continue; } lastW = w; lastQ = q; } return lastW; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :