--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowDepthCalculation.java	Tue Feb 27 18:06:52 2018 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowDepthCalculation.java	Wed Feb 28 17:27:15 2018 +0100
@@ -9,49 +9,33 @@
  */
 package org.dive4elements.river.artifacts.sinfo.flowdepth;
 
-import java.util.ArrayList;
-import java.util.Calendar;
 import java.util.Collection;
-import java.util.Collections;
-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.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.WKms;
 import org.dive4elements.river.artifacts.resources.Resources;
 import org.dive4elements.river.artifacts.sinfo.SINFOArtifact;
+import org.dive4elements.river.artifacts.sinfo.common.RiverInfoProvider;
 import org.dive4elements.river.artifacts.sinfo.flowdepth.FlowDepthAccess.DifferencesPair;
-import org.dive4elements.river.artifacts.sinfo.util.BedHeightInfo;
+import org.dive4elements.river.artifacts.sinfo.tkhcalculation.DischargeValuesFinder;
+import org.dive4elements.river.artifacts.sinfo.tkhcalculation.TkhCalculator;
+import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder;
+import org.dive4elements.river.artifacts.sinfo.util.CalculationUtils;
 import org.dive4elements.river.artifacts.sinfo.util.RiverInfo;
 import org.dive4elements.river.artifacts.sinfo.util.WstInfo;
 import org.dive4elements.river.artifacts.states.WaterlevelData;
 import org.dive4elements.river.artifacts.states.WaterlevelFetcher;
-import org.dive4elements.river.model.BedHeight;
-import org.dive4elements.river.model.BedHeightValue;
-import org.dive4elements.river.model.Gauge;
 import org.dive4elements.river.model.River;
-import org.dive4elements.river.utils.DoubleUtil;
-import org.dive4elements.river.utils.GaugeIndex;
 import org.dive4elements.river.utils.RiverUtils;
 
 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";
-
     private final CallContext context;
 
     public FlowDepthCalculation(final CallContext context) {
@@ -60,12 +44,7 @@
 
     public CalculationResult calculate(final SINFOArtifact sinfo) {
 
-        /*
-         * find the user of this artifact, sadly this is not part of the calling context, so instead we determine the
-         * owner oft the artifact
-         */
-        final ArtifactDatabase database = this.context.getDatabase();
-        final String user = database.findArtifactUser(sinfo.identifier());
+        final String user = CalculationUtils.findArtifactUser(this.context, sinfo);
 
         /* access input data */
         final FlowDepthAccess access = new FlowDepthAccess(sinfo);
@@ -74,24 +53,21 @@
 
         final Collection<DifferencesPair> diffPairs = access.getDifferencePairs();
 
-        final double from = access.getFrom();
-        final double to = access.getTo();
-        final DoubleRange calcRange = new DoubleRange(from, to);
+        final DoubleRange calcRange = access.getRange();
 
         final boolean useTkh = access.isUseTransportBodies();
 
         /* calculate results for each diff pair */
         final Calculation problems = new Calculation();
 
-        final List<Gauge> gauges = river.determineGauges(from, to);
-        final GaugeIndex gaugeIndex = new GaugeIndex(gauges);
+        final RiverInfoProvider infoProvider = RiverInfoProvider.forRange(this.context, river, calcRange);
 
         final String calcModeLabel = Resources.getMsg(this.context.getMeta(), sinfo.getCalculationMode().name());
 
         final FlowDepthCalculationResults results = new FlowDepthCalculationResults(calcModeLabel, user, riverInfo, calcRange, useTkh);
 
         for (final DifferencesPair diffPair : diffPairs) {
-            final FlowDepthCalculationResult result = calculateResult(river, calcRange, diffPair, problems, gaugeIndex, useTkh);
+            final FlowDepthCalculationResult result = calculateResult(calcRange, diffPair, problems, infoProvider, useTkh);
             if (result != null)
                 results.addResult(result);
         }
@@ -101,15 +77,17 @@
 
     /**
      * Calculates one W-MSH differences pair.
+     *
+     * @param infoProvider
      */
-    private FlowDepthCalculationResult calculateResult(final River river, final DoubleRange calcRange, final DifferencesPair diffPair,
-            final Calculation problems, final GaugeIndex gaugeIndex, final boolean useTkh) {
+    private FlowDepthCalculationResult calculateResult(final DoubleRange calcRange, final DifferencesPair diffPair,
+            final Calculation problems, final RiverInfoProvider infoProvider, final boolean useTkh) {
 
         /* access real input data from database */
         final String soundingId = diffPair.getSoundingId();
         final String wstId = diffPair.getWstId();
 
-        final BedHeight bedHeight = loadBedHeight(soundingId);
+        final BedHeightsFinder bedHeight = loadBedHeight(soundingId, calcRange);
         if (bedHeight == null) {
             final String message = Resources.format(this.context.getMeta(), "Failed to access sounding with id '{0}'", soundingId);
             problems.addProblem(message);
@@ -126,223 +104,29 @@
         final WKms wstKms = waterlevel.getWkms();
 
         final String wspLabel = wstKms.getName();
-        final String soundingLabel = bedHeight.getDescription();
+        final String soundingLabel = bedHeight.getInfo().getDescription();
         final String label = String.format("%s - %s", wspLabel, soundingLabel);
 
-        checkYearDifference(label, waterlevel, bedHeight, problems);
+        checkYearDifference(label, waterlevel, bedHeight.getInfo().getYear(), problems);
         checkWaterlevelDiscretisation(wstKms, calcRange, problems);
         // TODO: prüfen, ob sohlhöhen die calcRange abdecken/überschneiden
 
         /* re-determine the reference gauge, in the same way as the WaterlevelArtifact would do it */
-        final String notinrange = Resources.getMsg(this.context.getMeta(), CSV_NOT_IN_GAUGE_RANGE, CSV_NOT_IN_GAUGE_RANGE);
-
-        final Gauge refGauge = waterlevel.findReferenceGauge(river);
-        final String refGaugeName = refGauge == null ? notinrange : refGauge.getName();
-
-        final BedHeightInfo sounding = BedHeightInfo.from(bedHeight);
-        final int wspYear = waterlevel.getYear();
-        final WstInfo wstInfo = new WstInfo(wspLabel, wspYear, refGaugeName);
-
-        final FlowDepthCalculationResult resultData = new FlowDepthCalculationResult(label, wstInfo, sounding);
-
-        boolean doCalcTkh = useTkh;
-        if (doCalcTkh && !(wstKms instanceof QKms)) {
-            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
-            problems.addProblem(message);
-            doCalcTkh = false;
-        }
-
-        BedQualityD50KmValueFinder bedMeasurementsFinder = null;
-        if (doCalcTkh) {
-            bedMeasurementsFinder = loadBedMeasurements(river, calcRange, sounding.getYear().intValue());
-            if (bedMeasurementsFinder == null) {
-                final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
-                problems.addProblem(message);
-                doCalcTkh = false;
-            }
-        }
-
-        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 (doCalcTkh) {
-            qInterpolator = DoubleUtil.getLinearInterpolator(((QKms) wstKms).allKms(), ((QKms) wstKms).allQs());
-            if (qInterpolator != null)
-                qRange = new DoubleRange(((QKms) wstKms).allQs().min(), ((QKms) wstKms).allQs().max());
-            else {
-                final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
-                problems.addProblem(message);
-                doCalcTkh = false;
-            }
-        }
-
-        // 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());
-
-        // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
-        /* SoilKind lastKind = SoilKind.mobil; */
-        SoilKindKmValueFinder soilKindFinder = null;
-        if (doCalcTkh) {
-            soilKindFinder = new SoilKindKmValueFinder();
-            if (!soilKindFinder.loadValues(river, calcRange)) {
-                doCalcTkh = false;
-                final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
-                problems.addProblem(message);
-            }
-        }
-
-        FlowVelocityModelKmValueFinder flowVelocitiesFinder = null;
-        if (doCalcTkh) {
-            flowVelocitiesFinder = new FlowVelocityModelKmValueFinder();
-            if (!flowVelocitiesFinder.loadValues(river, calcRange, qRange)) {
-                doCalcTkh = false;
-                final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, label);
-                problems.addProblem(message);
-            }
-        }
-
-        for (final BedHeightValue bedHeightValue : sortedValues) {
-
-            final Double station = bedHeightValue.getStation();
-            if (station == null || station.isNaN())
-                continue;
-
-            final Double meanBedHeightDbl = bedHeightValue.getHeight();
-            if (meanBedHeightDbl == null || meanBedHeightDbl.isNaN())
-                continue;
-
-            final double km = station;
-            final double meanBedHeight = meanBedHeightDbl;
-
-            if (!calcRange.containsDouble(km))
-                continue;
+        final RiverInfoProvider riverInfoProvider = infoProvider.forWaterlevel(waterlevel);
 
-            try {
-                // FIXME: check out of range
-                final double wst = wstInterpolator.value(km);
-
-                final double flowDepth = wst - meanBedHeight;
-
-                // FIXME: piecewise constant interpolation?
-                // final double discharge = wstKms instanceof QKms ? ((QKms) wstKms).getQ(i) : Double.NaN;
-                double discharge;
-                if (qInterpolator != null)
-                    discharge = qInterpolator.value(km);
-                else
-                    discharge = Double.NaN;
-
-                // Calculate tkh
-                double tkh = 0;
-                if (doCalcTkh) {
-                    double d50 = Double.NaN;
-                    try {
-                        d50 = bedMeasurementsFinder.findD50(km);
-                    }
-                    catch (final 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 to one message?
-                    }
-                    if (!Double.isNaN(d50)) {
-                        if (flowVelocitiesFinder.findKmQValues(km, discharge)) {
-                            tkh = calculateTkh(wst - meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), d50, flowVelocitiesFinder.getFindTauFound());
-                            if (!Double.isNaN(tkh) && (tkh < 0))
-                                tkh = 0;
-                            /*
-                             * 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 {
-                            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
-                            problems.addProblem(km, message);
-                            // FIXME: cumulate problems to one message?
-                        }
-                    } else
-                        tkh = Double.NaN;
-                }
+        final int wspYear = waterlevel.getYear();
+        final WstInfo wstInfo = new WstInfo(wspLabel, wspYear, riverInfoProvider.getReferenceGauge());
 
-                // Soil kind
-                SoilKind kind = SoilKind.mobil;
-                if (doCalcTkh) {
-                    try {
-                        kind = soilKindFinder.findSoilKind(km);
-                    }
-                    catch (final 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 to one message?
-                    }
-                }
-
-                final double flowDepthTkh;
-                final double tkhUp;
-                final double tkhDown;
-                switch (kind) {
-                case starr:
-                    flowDepthTkh = wst - (meanBedHeight + tkh / 100);
-                    tkhUp = tkh;
-                    tkhDown = 0;
-                    break;
+        final DischargeValuesFinder dischargeProvider = DischargeValuesFinder.fromKms(wstKms);
 
-                case mobil:
-                default:
-                    flowDepthTkh = wst - (meanBedHeight + tkh / 200);
-                    tkhUp = tkh / 2;
-                    tkhDown = -tkh / 2;
-                    break;
-                }
-
-                // REMARK: access the location once only during calculation
-                final String location = LocationProvider.getLocation(river.getName(), km);
+        final River river = riverInfoProvider.getRiver();
+        final TkhCalculator tkhCalculator = TkhCalculator.buildTkhCalculator(useTkh, this.context, problems, label, river, calcRange, dischargeProvider,
+                bedHeight);
 
-                // REMARK: access the gauge once only during calculation
-                final Gauge gauge = findGauge(waterlevel, refGauge, gaugeIndex, km);
-
-                final String gaugeLabel = gauge == null ? notinrange : gauge.getName();
-
-                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 */
-                e.printStackTrace();
-                /* simply ignore */
-            }
-        }
-
-        return resultData;
+        final FlowDepthCalculator calculator = new FlowDepthCalculator(riverInfoProvider, wstKms, dischargeProvider, bedHeight, tkhCalculator);
+        return calculator.execute(label, wstInfo, calcRange);
     }
 
-    /**
-     * Sohlbeschaffenheit (D50 Korndurchmesser aus Seddb)
-     * Abhängig von Peiljahr
-     */
-    private BedQualityD50KmValueFinder loadBedMeasurements(final River river, final DoubleRange kmRange, final int soundingYear) {
-
-        /* construct valid measurement time range */
-        final Calendar cal = Calendar.getInstance();
-        cal.clear();
-
-        cal.set(soundingYear - VALID_BED_MEASUREMENT_YEARS, 0, 1);
-        final Date startTime = cal.getTime();
-
-        cal.set(soundingYear + VALID_BED_MEASUREMENT_YEARS, 11, 31);
-        final Date endTime = cal.getTime();
-
-        final BedQualityD50KmValueFinder finder = new BedQualityD50KmValueFinder();
-        if (finder.loadValues(river, kmRange, new DateRange(startTime, endTime)))
-            return finder;
-        else
-            return null;
-    }
 
     /**
      * Checks the year difference between waterlevels and sounding, and issues a warning if too big.
@@ -353,9 +137,7 @@
      * 1918 ≤ X < 1958 ± 12
      * X < 1918 ± 25
      */
-    private void checkYearDifference(final String label, final WaterlevelData waterlevel, final BedHeight sounding, final Calculation problems) {
-
-        final Integer soundingYear = sounding.getYear();
+    private void checkYearDifference(final String label, final WaterlevelData waterlevel, final Integer soundingYear, final Calculation problems) {
         if (soundingYear == null)
             return;
 
@@ -388,21 +170,6 @@
         return 3;
     }
 
-    private Gauge findGauge(final WaterlevelData waterlevel, final Gauge refGauge, final GaugeIndex gaugeIndex, final double km) {
-
-        // REMARK: using same logic as in WaterlevelExporter here
-
-        final boolean showAllGauges = waterlevel.isShowAllGauges();
-
-        if (showAllGauges)
-            return gaugeIndex.findGauge(km);
-
-        if (refGauge.getRange().contains(km))
-            return refGauge;
-
-        return null;
-    }
-
     /* Checks if the discretisation of the waterlevel exceeds 1000m */
 
     private void checkWaterlevelDiscretisation(final WKms wstKms, final DoubleRange calcRange, final Calculation problems) {
@@ -424,7 +191,7 @@
         }
     }
 
-    private BedHeight loadBedHeight(final String soundingId) {
+    private BedHeightsFinder loadBedHeight(final String soundingId, final DoubleRange calcRange) {
 
         // REMARK: absolutely unbelievable....
         // The way how bed-heights (and other data too) is accessed is different for nearly every calculation-type
@@ -438,6 +205,11 @@
         final BedHeightsArtifact artifact = (BedHeightsArtifact) RiverUtils.getArtifact(parts[0], this.context);
 
         final Integer bedheightId = artifact.getDataAsInteger("height_id");
+        if (bedheightId == null) {
+            // FIXME: error message!
+            return null;
+        }
+
         // REMARK: this only works with type 'single'; unclear on how to distinguish from epoch data (or whatever the
         // other type means)
         // Luckily, the requirement is to only access 'single' data here.
@@ -449,35 +221,6 @@
         // BedHeightFactory uses its own (direct) way of accessing the data, with its own implemented data classes.
         // return BedHeightFactory.getHeight(bedheightType, bedheightId, from, to);
 
-        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(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);
+        return BedHeightsFinder.forId(bedheightId, calcRange);
     }
 }
\ No newline at end of file