/** 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.dive4elements.river.artifacts.model.Calculation;
import org.dive4elements.river.artifacts.sinfo.common.SInfoResultRow;
import org.dive4elements.river.artifacts.sinfo.common.SInfoResultType;
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 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 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(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);

        if (!dischargeProvider.isValid()) {
            problems.addProblem("sinfo_calc_flow_depth.warning.missingQ", label);
            return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
        }

        /* access bed quality data */
        final int soundingYear = bedHeightsProvider.getInfo().getYear();
        final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(problems, river, calcRange, soundingYear,
                VALID_BED_MEASUREMENT_YEARS);
        if (bedMeasurementsFinder == null)
            return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);

        /* access bed soil kind data */
        final SoilKindKmValueFinder soilKindFinder = SoilKindKmValueFinder.loadValues(problems, river, calcRange);
        if (soilKindFinder == null)
            return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);

        final DoubleRange qRange = dischargeProvider.getRange();
        final FlowVelocityModelKmValueFinder flowVelocitiesFinder = FlowVelocityModelKmValueFinder.loadValues(problems, river, calcRange, qRange);
        if (flowVelocitiesFinder == null)
            return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);

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

    private TkhCalculator(final BedQualityD50KmValueFinder bedMeasurementsFinder, final WaterlevelValuesFinder waterlevelProvider,
            final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider, final SoilKindKmValueFinder soilKindFinder,
            final FlowVelocityModelKmValueFinder flowVelocitiesFinder) {
        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) {

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

        return this.soilKindFinder.findSoilKind(km);
    }

    private double getBedMeasurement(final double km) {
        return this.bedMeasurementsFinder.findD50(km);
    }

    public boolean calculateTkh(final double km, final SInfoResultRow row) {

        row.putValue(SInfoResultType.station, km);

        final SoilKind kind = getSoilKind(km);
        row.putValue(SInfoResultType.soilkind, kind);

        final double wst = this.waterlevelProvider.getWaterlevel(km);
        row.putValue(SInfoResultType.waterlevel, wst);
        if (Double.isNaN(wst))
            return false;

        final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km);
        row.putValue(SInfoResultType.meanBedHeight, meanBedHeight);
        if (Double.isNaN(meanBedHeight))
            return false;

        final double flowDepth = wst - meanBedHeight;
        row.putValue(SInfoResultType.flowdepth, flowDepth);
        if (Double.isNaN(flowDepth))
            return false;

        final double discharge = this.dischargeProvider.getDischarge(km);
        row.putValue(SInfoResultType.discharge, discharge);

        if (!this.hasTkh())
            return true;

        // Missing discharge or kind is only a problem if we want to calculate tkh
        if (Double.isNaN(discharge))
            return false;
        if (kind == null)
            return false;

        final double d50 = getBedMeasurement(km);
        row.putValue(SInfoResultType.d50, d50);
        if (Double.isNaN(d50))
            return false;

        if (!this.flowVelocitiesFinder.findKmQValues(km, discharge))
            return false;

        final double velocity = this.flowVelocitiesFinder.getFindVmainFound();
        row.putValue(SInfoResultType.velocity, velocity);
        if (Double.isNaN(velocity))
            return false;

        final double tau = this.flowVelocitiesFinder.getFindTauFound();
        row.putValue(SInfoResultType.tau, tau);
        if (Double.isNaN(tau))
            return false;

        final double tkh = calculateTkh(wst - meanBedHeight, velocity, d50, tau);
        row.putValue(SInfoResultType.tkh, tkh);
        if (Double.isNaN(tkh))
            return false;

        switch (kind) {
        case starr:
            row.putValue(SInfoResultType.tkhup, tkh);
            row.putValue(SInfoResultType.tkhdown, 0.0);
            break;

        case mobil:
        default:
            row.putValue(SInfoResultType.tkhup, tkh / 2);
            row.putValue(SInfoResultType.tkhdown, -tkh / 2);
            break;
        }

        final double flowDepthTkh = calculateFlowDepthTkh(tkh, kind, wst, meanBedHeight);
        row.putValue(SInfoResultType.flowdepthtkh, flowDepthTkh);

        return true;
    }

    /**
     * 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) {

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

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