/** 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.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 WaterlevelValuesFinder waterlevelProvider;

    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 WaterlevelValuesFinder waterlevelProvider, final DischargeValuesFinder dischargeProvider,
            final BedHeightsFinder bedHeightsProvider) {

        if (!useTkh)
            return new TkhCalculator(problems, label, context, null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);

        if (!dischargeProvider.isValid()) {
            final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
            problems.addProblem(message);
            return new TkhCalculator(problems, label, context, null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
        }

        final int 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.missingD50", null, label);
            problems.addProblem(message);
            return new TkhCalculator(problems, label, context, null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, 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 new TkhCalculator(problems, label, context, null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, 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 new TkhCalculator(problems, label, context, null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
        }

        return new TkhCalculator(problems, label, context, bedMeasurementsFinder, waterlevelProvider, dischargeProvider, bedHeightsProvider, soilKindFinder,
                flowVelocitiesFinder);
    }

    private TkhCalculator(final Calculation problems, final String problemLabel, final CallContext context,
            final BedQualityD50KmValueFinder bedMeasurementsFinder, final WaterlevelValuesFinder waterlevelProvider,
            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.waterlevelProvider = waterlevelProvider;
        this.dischargeProvider = dischargeProvider;
        this.bedHeightsProvider = bedHeightsProvider;
        this.soilKindFinder = soilKindFinder;
        this.flowVelocitiesFinder = flowVelocitiesFinder;
    }

    public boolean hasTkh() {

        if (this.dischargeProvider == null || !this.dischargeProvider.isValid())
            return false;

        if (this.bedMeasurementsFinder == null)
            return false;

        if (this.soilKindFinder == null)
            return false;

        if (this.flowVelocitiesFinder == null)
            return false;

        return true;
    }

    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 SoilKind kind = getSoilKind(km);

        final double wst = this.waterlevelProvider.getWaterlevel(km);

        final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km);

        final double flowDepth = wst - meanBedHeight;

        final double discharge = this.dischargeProvider.getDischarge(km);
        if (Double.isNaN(discharge))
            return new Tkh(km, wst, meanBedHeight, flowDepth, Double.NaN, kind);

        if (!this.hasTkh())
            return new Tkh(km, wst, meanBedHeight, flowDepth, Double.NaN, kind);

        final double d50 = getBedMeasurement(km);
        if (Double.isNaN(d50))
            return new Tkh(km, wst, meanBedHeight, flowDepth, discharge, kind);

        if (!this.flowVelocitiesFinder.findKmQValues(km, discharge)) {
            // TODO: ggf. station in Fehlermeldung?
            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, this.problemLabel);
            this.problems.addProblem(km, message);
            // FIXME: cumulate problems to one message?
            return new Tkh(km, wst, meanBedHeight, flowDepth, discharge, kind);
        }

        final double velocity = this.flowVelocitiesFinder.getFindVmainFound();
        final double tau = this.flowVelocitiesFinder.getFindTauFound();

        final double tkh = calculateTkh(wst - meanBedHeight, velocity, d50, tau);
        double tkhUp;
        double tkhDown;
        switch (kind) {
        case starr:
            tkhUp = tkh;
            tkhDown = 0;
            break;

        case mobil:
        default:
            tkhUp = tkh / 2;
            tkhDown = -tkh / 2;
            break;
        }

        final double flowDepthTkh = calculateFlowDepthTkh(tkhUp, kind, wst, meanBedHeight);

        return new Tkh(km, wst, meanBedHeight, flowDepth, flowDepthTkh, discharge, kind, tkh, tkhUp, tkhDown, velocity, d50, tau);
    }

    /**
     * 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 (!), never negative
     */
    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;
        final double tkh = 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
        // Some regular input values may give a negative calculation result; that is unwanted
        if (tkh < 0.0)
            return 0.0;

        return tkh;
    }

    private double calculateFlowDepthTkh(final double tkhValue, final SoilKind tkhKind, final double wst, final double meanBedHeight) {

        if (Double.isNaN(tkhValue) || tkhKind == null)
            return Double.NaN;

        switch (tkhKind) {
        case starr:
            return wst - (meanBedHeight + tkhValue / 100);

        case mobil:
        default:
            return wst - (meanBedHeight + tkhValue / 200);
        }
    }
}