--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/TkhCalculator.java	Wed Feb 28 17:27:15 2018 +0100
@@ -0,0 +1,232 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ *  Björnsen Beratende Ingenieure GmbH
+ *  Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.sinfo.tkhcalculation;
+
+import org.apache.commons.lang.math.DoubleRange;
+import org.apache.commons.math.ArgumentOutsideDomainException;
+import org.apache.commons.math.FunctionEvaluationException;
+import org.dive4elements.artifacts.CallContext;
+import org.dive4elements.river.artifacts.model.Calculation;
+import org.dive4elements.river.artifacts.resources.Resources;
+import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder;
+import org.dive4elements.river.model.River;
+
+/**
+ * @author Gernot Belger
+ */
+public final class TkhCalculator {
+
+    private static final int VALID_BED_MEASUREMENT_YEARS = 20;
+
+    private final Calculation problems;
+
+    private final String problemLabel;
+
+    private final CallContext context;
+
+    private final BedQualityD50KmValueFinder bedMeasurementsFinder;
+
+    private final SoilKindKmValueFinder soilKindFinder;
+
+    private final BedHeightsFinder bedHeightsProvider;
+
+    private final DischargeValuesFinder dischargeProvider;
+
+    private final FlowVelocityModelKmValueFinder flowVelocitiesFinder;
+
+    public static TkhCalculator buildTkhCalculator(final boolean useTkh, final CallContext context, final Calculation problems, final String label,
+            final River river, final DoubleRange calcRange, final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider) {
+
+        if (!useTkh)
+            return null;
+
+        if (!dischargeProvider.isValid()) {
+            final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
+            problems.addProblem(message);
+            return null;
+        }
+
+        final Integer soundingYear = bedHeightsProvider.getInfo().getYear();
+        final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(river, calcRange, soundingYear,
+                VALID_BED_MEASUREMENT_YEARS);
+
+        if (bedMeasurementsFinder == null) {
+            final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
+            problems.addProblem(message);
+            return null;
+        }
+
+        // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
+        final SoilKindKmValueFinder soilKindFinder = SoilKindKmValueFinder.loadValues(river, calcRange);
+        if (soilKindFinder == null) {
+            final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
+            problems.addProblem(message);
+            return null;
+        }
+
+        final DoubleRange qRange = dischargeProvider.getRange();
+        final FlowVelocityModelKmValueFinder flowVelocitiesFinder = FlowVelocityModelKmValueFinder.loadValues(river, calcRange, qRange);
+        if (flowVelocitiesFinder == null) {
+            final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, label);
+            problems.addProblem(message);
+            return null;
+        }
+
+        return new TkhCalculator(problems, label, context, bedMeasurementsFinder, dischargeProvider, bedHeightsProvider, soilKindFinder, flowVelocitiesFinder);
+    }
+
+    private TkhCalculator(final Calculation problems, final String problemLabel, final CallContext context,
+            final BedQualityD50KmValueFinder bedMeasurementsFinder, final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider,
+            final SoilKindKmValueFinder soilKindFinder,
+            final FlowVelocityModelKmValueFinder flowVelocitiesFinder) {
+        this.problems = problems;
+        this.problemLabel = problemLabel;
+        this.context = context;
+        this.bedMeasurementsFinder = bedMeasurementsFinder;
+        this.dischargeProvider = dischargeProvider;
+        this.bedHeightsProvider = bedHeightsProvider;
+        this.soilKindFinder = soilKindFinder;
+        this.flowVelocitiesFinder = flowVelocitiesFinder;
+    }
+
+    private double getDischarge(final double km) {
+
+        try {
+            return this.dischargeProvider.getDischarge(km);
+        }
+        catch (final FunctionEvaluationException e) {
+            // TODO: exceptions nicht komplett schlucken? evtl. mit log.debug(e) ausgeben
+            return Double.NaN;
+        }
+    }
+
+    private SoilKind getSoilKind(final double km) {
+
+        try {
+            return this.soilKindFinder.findSoilKind(km);
+        }
+        catch (final ArgumentOutsideDomainException e) {
+            // FIXME: cumulate problems to one message?
+            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, this.problemLabel);
+            this.problems.addProblem(km, message);
+            return null;
+        }
+    }
+
+    private double getBedMeasurement(final double km) {
+
+        try {
+            return this.bedMeasurementsFinder.findD50(km);
+        }
+        catch (final Exception e) {
+            // FIXME: cumulate problems to one message?
+            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, this.problemLabel);
+            this.problems.addProblem(km, message);
+
+            return Double.NaN;
+        }
+    }
+
+    public Tkh getTkh(final double km, final double wst) {
+
+        final SoilKind kind = getSoilKind(km);
+
+        final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km);
+
+        final double discharge = getDischarge(km);
+        if (Double.isNaN(discharge)) {
+
+            // final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null,
+            // this.problemLabel);
+            // this.problems.addProblem(km, message);
+
+            // TODO: nochmal gemeinsam überlegen welche probleme wir loggen, an dieser stelle müsste man ggf. die station
+            // mitausgeben
+
+            return new Tkh(km, wst, meanBedHeight, Double.NaN, kind, Double.NaN, Double.NaN, Double.NaN);
+        }
+
+        final double d50 = getBedMeasurement(km);
+        if (Double.isNaN(d50))
+            return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
+
+        if (!this.flowVelocitiesFinder.findKmQValues(km, discharge)) {
+            // TODO: ggf. station in Fehlermeldung?
+            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel);
+            this.problems.addProblem(km, message);
+            // FIXME: cumulate problems to one message?
+        }
+
+        final double tkh = calculateTkh(wst - meanBedHeight, this.flowVelocitiesFinder.getFindVmainFound(), d50, this.flowVelocitiesFinder.getFindTauFound());
+        // FIXME: noch mal prüfen, im alten code wurde hier immer auf 0 gesetzt
+        if (Double.isNaN(tkh) || (tkh < 0)) {
+            // TODO: ggf. station in Fehlermeldung?
+
+            // FIXME: Fehlermeldung nicht korrekt, passiert mit Wasserspiegel 'MHQ' und 'QP-1993': alle Daten (auch Abfluss)
+            // vorhanden, aber tkh negativ...
+            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel);
+            this.problems.addProblem(km, message);
+
+            return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
+        }
+
+        /*
+         * log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f",
+         * km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(),
+         * d50*1000, tkh));
+         */
+
+        double tkhUp;
+        double tkhDown;
+        switch (kind) {
+        case starr:
+            tkhUp = tkh;
+            tkhDown = 0;
+            break;
+
+        case mobil:
+        default:
+            tkhUp = tkh / 2;
+            tkhDown = -tkh / 2;
+            break;
+        }
+
+        return new Tkh(km, wst, meanBedHeight, discharge, kind, tkh, tkhUp, tkhDown);
+    }
+
+    /**
+     * Calculates a transport body height
+     *
+     * @param h
+     *            flow depth in m
+     * @param vm
+     *            flow velocity in m
+     * @param d50
+     *            grain diameter D50 in m (!)
+     * @param tau
+     *            shear stress in N/m^2
+     * @return transport body height in cm (!)
+     */
+    private double calculateTkh(final double h, final double vm, final double d50, final double tau) {
+        final double PHYS_G = 9.81;
+        final double PHYS_SPECGRAV_S = 2.6;
+        final double PHYS_VELOCCOEFF_N = 6;
+        final double PHYS_FORMCOEFF_ALPHA = 0.7;
+        final double PHYS_VISCOSITY_NUE = 1.3e-6;
+        final double PHYS_GRAIN_DENSITY_RHOS = 2603;
+        final double PHYS_WATER_DENSITY_RHO = 999.97;
+
+        final double froude = vm / Math.sqrt(PHYS_G * h);
+        final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50;
+        final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6));
+        final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50;
+        return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
+    }
+}
\ No newline at end of file