/** 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);
    }
}