--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowDepthCalculation.java	Wed Feb 14 19:06:21 2018 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowDepthCalculation.java	Thu Feb 22 12:03:31 2018 +0100
@@ -16,18 +16,19 @@
 import java.util.Date;
 import java.util.List;
 
+import org.apache.commons.lang.math.DoubleRange;
 import org.apache.commons.math.FunctionEvaluationException;
 import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.log4j.Logger;
 import org.dive4elements.artifacts.ArtifactDatabase;
 import org.dive4elements.artifacts.CallContext;
 import org.dive4elements.river.artifacts.BedHeightsArtifact;
 import org.dive4elements.river.artifacts.model.Calculation;
 import org.dive4elements.river.artifacts.model.CalculationResult;
+import org.dive4elements.river.artifacts.model.DateRange;
 import org.dive4elements.river.artifacts.model.LocationProvider;
-import org.dive4elements.river.artifacts.model.QKms;
+import org.dive4elements.river.artifacts.model.WQKms;
 import org.dive4elements.river.artifacts.model.WKms;
-import org.dive4elements.river.artifacts.model.minfo.QualityMeasurementFactory;
-import org.dive4elements.river.artifacts.model.minfo.QualityMeasurements;
 import org.dive4elements.river.artifacts.resources.Resources;
 import org.dive4elements.river.artifacts.sinfo.SINFOArtifact;
 import org.dive4elements.river.artifacts.sinfo.flowdepth.FlowDepthAccess.DifferencesPair;
@@ -43,6 +44,8 @@
 
 class FlowDepthCalculation {
 
+    private static Logger log = Logger.getLogger(FlowDepthCalculation.class);
+
     private static final int VALID_BED_MEASUREMENT_YEARS = 20;
 
     private static final String CSV_NOT_IN_GAUGE_RANGE = "export.waterlevel.csv.not.in.gauge.range";
@@ -70,6 +73,7 @@
 
         final double from = access.getFrom();
         final double to = access.getTo();
+        final DoubleRange calcRange = new DoubleRange(from, to);
 
         final boolean useTkh = access.isUseTransportBodies();
 
@@ -84,7 +88,7 @@
         final FlowDepthCalculationResults results = new FlowDepthCalculationResults(calcModeLabel, user, river, from, to, useTkh);
 
         for (final DifferencesPair diffPair : diffPairs) {
-            final FlowDepthCalculationResult result = calculateResult(river, from, to, diffPair, problems, gaugeIndex);
+            final FlowDepthCalculationResult result = calculateResult(river, calcRange, diffPair, problems, gaugeIndex, useTkh);
             if (result != null)
                 results.addResult(result);
         }
@@ -92,14 +96,17 @@
         return new CalculationResult(results, problems);
     }
 
-    private FlowDepthCalculationResult calculateResult(final River river, final double from, final double to, final DifferencesPair diffPair,
-            final Calculation problems, final GaugeIndex gaugeIndex) {
+    /**
+     * Calculates one W-MSH differences pair.
+     */
+    private FlowDepthCalculationResult calculateResult(final River river, final DoubleRange calcRange, final DifferencesPair diffPair,
+            final Calculation problems, final GaugeIndex gaugeIndex, final boolean useTkh) {
 
         /* access real input data from database */
         final String soundingId = diffPair.getSoundingId();
         final String wstId = diffPair.getWstId();
 
-        final BedHeight bedHeight = loadBedHeight(soundingId, from, to);
+        final BedHeight bedHeight = loadBedHeight(soundingId);
         if (bedHeight == null) {
             final String message = Resources.format(this.context.getMeta(), "Failed to access sounding with id '{0}'", soundingId);
             problems.addProblem(message);
@@ -136,36 +143,49 @@
 
         // FIXME: nur prüfen/beschaffen wenn TKH Berechnung aktiv
         /* Abflusswerte vorhanden? */
-        if (!(wstKms instanceof QKms)) {
+        if (!(wstKms instanceof WQKms)) {
             final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
             problems.addProblem(message);
             // TODO: keine Berechnung TKH
         }
 
-        final QualityMeasurements bedMeasurements = getBedMeasurements(river, from, to, sounding.getYear());
-        // FIXME: prüfung ob (genug) werte vorhanden sind? was sind genau die kriterien? falls nein, problemhinzufügen und keine
+        BedQualityD50KmValueFinder bedMeasurementsFinder = null;
+        if (useTkh)
+            bedMeasurementsFinder = loadBedMeasurements(river, calcRange, sounding.getYear());
+        // FIXME: prüfung ob (genug) werte vorhanden sind? was sind genau die kriterien? falls nein, problem hinzufügen und keine
         // berechnung tkh
-        // FIXME: wie wird ggf. interpoliert? --> absprache?
-        // FIXME: mir direkt aufgefallen, die Beispieldatenbank liefert Werte zum gleichen km und zeitpunkt, die messwerte sind
-        // aber unterschiedlich....???
-        // FIXME: die eigentlichen daten extrahieren, ggf. wenn esswerte zum gleichen datum voriliegen. das neueste nehmen? oder
-        // das näheste zum Peiljahr?
-
-        // FIXME Art der Gewässersohle (starr/mobil)
-        // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
 
         final String bedHeightLabel = bedHeight.getDescription();
         final String wstLabel = wstKms.getName();
 
         final UnivariateRealFunction wstInterpolator = DoubleUtil.getLinearInterpolator(wstKms.allKms(), wstKms.allWs());
-
+        UnivariateRealFunction qInterpolator = null;
+        DoubleRange qRange = null;
+        if (useTkh && (wstKms instanceof WQKms)) {
+            qInterpolator = DoubleUtil.getLinearInterpolator(((WQKms) wstKms).allKms(), ((WQKms) wstKms).allQs());
+            qRange = new DoubleRange( ((WQKms) wstKms).allQs().min(), ((WQKms) wstKms).allQs().max());
+        }
+        
         // FIXME: sort by station first, but in what direction?
+        // FIXME: using river.getKmUp()?
         final List<BedHeightValue> values = bedHeight.getValues();
 
         final List<BedHeightValue> sortedValues = new ArrayList<>(values);
         Collections.sort(sortedValues, new BedHeightStationComparator());
 
-        SoilKind lastKind = SoilKind.mobil;
+        // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
+        /* SoilKind lastKind = SoilKind.mobil; */
+        SoilKindKmValueFinder soilKindFinder = null;
+        if (useTkh) {
+            soilKindFinder = new SoilKindKmValueFinder();
+            soilKindFinder.loadValues(river, calcRange);
+        }
+        
+        FlowVelocityModelKmValueFinder flowVelocitiesFinder = null;
+        if (useTkh) {
+            flowVelocitiesFinder = new FlowVelocityModelKmValueFinder();
+            flowVelocitiesFinder.loadValues(river, calcRange, qRange);
+        }
 
         for (final BedHeightValue bedHeightValue : sortedValues) {
 
@@ -180,6 +200,9 @@
             final double km = station;
             final double meanBedHeight = meanBedHeightDbl;
 
+            if (!calcRange.containsDouble(km))
+                continue;
+            
             try {
                 // FIXME: check out of range
                 final double wst = wstInterpolator.value(km);
@@ -188,30 +211,63 @@
 
                 // FIXME: piecewise constant interpolation?
                 // final double discharge = wstKms instanceof QKms ? ((QKms) wstKms).getQ(i) : Double.NaN;
-                final double discharge = Double.NaN;
+                double discharge;
+                if (qInterpolator != null)
+                    discharge = qInterpolator.value(km);
+                else
+                    discharge = Double.NaN;
 
                 // FIXME: calculate tkh
-
-                // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt
+                double tkh = 0;
+                if (useTkh) {
+                    double d50 = 0;
+                    try {
+                        d50 = bedMeasurementsFinder.findD50(km);
+                    } catch (Exception e) {
+                        final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, label);
+                        problems.addProblem(km, message);
+                        //FIXME: cumulate problems
+                    }
+                    if (flowVelocitiesFinder.findKmQValues(km, discharge)) {
+                        tkh = calculateTkh(wst - meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), d50, flowVelocitiesFinder.getFindTauFound());
+                        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));
+                    }
+                    else
+                        tkh = Double.NaN;
+                }
+                
+                // Soil kind
+                SoilKind kind = SoilKind.starr;
+                if (useTkh) {
+                    try {
+                        kind = soilKindFinder.findSoilKind(km);                        
+                    } catch (Exception e) {
+                        final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
+                        problems.addProblem(km, message);
+                        //FIXME: cumulate problems
+                    }
+                }
+                
+                /* // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt
+                SoilKind kind;
                 final boolean changeKind = Math.random() > 0.95;
-                SoilKind kind;
                 if (changeKind) {
                     switch (lastKind) {
                     case starr:
                         kind = SoilKind.mobil;
                         break;
-
                     case mobil:
                     default:
                         kind = SoilKind.starr;
                         break;
-
                     }
                 } else
                     kind = lastKind;
                 lastKind = kind;
-
-                final double tkh = 100 + 10 * (Math.random() - 0.5);
+                */
+                
+                //tkh = 100 + 10 * (Math.random() - 0.5);
 
                 final double flowDepthTkh;
                 final double tkhUp;
@@ -240,8 +296,8 @@
 
                 final String gaugeLabel = gauge == null ? notinrange : gauge.getName();
 
-                resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel, gaugeLabel, meanBedHeight, bedHeightLabel,
-                        location);
+                resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel,
+                        gaugeLabel, meanBedHeight, bedHeightLabel, location);
             }
             catch (final FunctionEvaluationException e) {
                 /* should only happen if out of range */
@@ -258,7 +314,7 @@
      * Sohlbeschaffenheit (D50 Korndurchmesser aus Seddb)
      * Abhängig von Peiljahr
      */
-    private QualityMeasurements getBedMeasurements(final River river, final double from, final double to, final int soundingYear) {
+    private BedQualityD50KmValueFinder loadBedMeasurements(final River river, final DoubleRange kmRange, final int soundingYear) {
 
         /* construct valid measurement time range */
         final Calendar cal = Calendar.getInstance();
@@ -270,7 +326,11 @@
         cal.set(soundingYear + VALID_BED_MEASUREMENT_YEARS, 11, 31);
         final Date endTime = cal.getTime();
 
-        return QualityMeasurementFactory.getBedMeasurements(river.getName(), from, to, startTime, endTime);
+        final BedQualityD50KmValueFinder finder = new BedQualityD50KmValueFinder();
+        if (finder.loadValues(river, kmRange, new DateRange(startTime, endTime)))
+            return finder;
+        else
+            return null;
     }
 
     /**
@@ -348,10 +408,10 @@
         }
     }
 
-    private BedHeight loadBedHeight(final String soundingId, final double from, final double to) {
+    private BedHeight loadBedHeight(final String soundingId) {
 
         // REMARK: absolutely unbelievable....
-        // The way how bed-heights (and other data too) is accessed is different for nearly ever calculation-type
+        // The way how bed-heights (and other data too) is accessed is different for nearly every calculation-type
         // throughout flys.
         // The knowledge on how to parse the datacage-ids is spread through the complete code-base...
 
@@ -375,4 +435,28 @@
 
         return BedHeight.getBedHeightById(bedheightId);
     }
+    
+    /**
+     * 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(double h, double vm, double d50, 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