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