/* 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.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.apache.log4j.Logger;
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.flowdepth.FlowDepthAccess.DifferencesPair;
import org.dive4elements.river.artifacts.sinfo.util.BedHeightInfo;
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) {
        this.context = context;
    }

    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());

        /* access input data */
        final FlowDepthAccess access = new FlowDepthAccess(sinfo);
        final River river = access.getRiver();
        final RiverInfo riverInfo = new RiverInfo(river);

        final Collection<DifferencesPair> diffPairs = access.getDifferencePairs();

        final double from = access.getFrom();
        final double to = access.getTo();
        final DoubleRange calcRange = new DoubleRange(from, to);

        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 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);
            if (result != null)
                results.addResult(result);
        }

        return new CalculationResult(results, problems);
    }

    /**
     * Calculates one W-MSH differences pair.
     */
    private FlowDepthCalculationResult calculateResult(final River river, final DoubleRange calcRange, final DifferencesPair diffPair,
            final Calculation problems, final GaugeIndex gaugeIndex, final boolean useTkh) {

        /* access real input data from database */
        final String soundingId = diffPair.getSoundingId();
        final String wstId = diffPair.getWstId();

        final BedHeight bedHeight = loadBedHeight(soundingId);
        if (bedHeight == null) {
            final String message = Resources.format(this.context.getMeta(), "Failed to access sounding with id '{0}'", soundingId);
            problems.addProblem(message);
            return null;
        }

        /* REMARK: fetch ALL wst kms, because we want to determine the original reference gauge */
        final WaterlevelData waterlevel = new WaterlevelFetcher().findWaterlevel(this.context, wstId, Double.NaN, Double.NaN);
        if (waterlevel == null) {
            final String message = Resources.format(this.context.getMeta(), "Failed to access waterlevel with id '{0}'", wstId);
            problems.addProblem(message);
            return null;
        }
        final WKms wstKms = waterlevel.getWkms();

        final String wspLabel = wstKms.getName();
        final String soundingLabel = bedHeight.getDescription();
        final String label = String.format("%s - %s", wspLabel, soundingLabel);

        checkYearDifference(label, waterlevel, bedHeight, 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;

            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 (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;
                }
                
                // Soil kind
                SoilKind kind = SoilKind.mobil;
                if (doCalcTkh) {
                    try {
                        kind = soilKindFinder.findSoilKind(km);                        
                    } catch (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?
                    }
                }
                
                /* // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt
                SoilKind kind;
                final boolean changeKind = Math.random() > 0.95;
                if (changeKind) {
                    switch (lastKind) {
                    case starr:
                        kind = SoilKind.mobil;
                        break;
                    case mobil:
                    default:
                        kind = SoilKind.starr;
                        break;
                    }
                } else
                    kind = lastKind;
                lastKind = kind;
                */
                
                //tkh = 100 + 10 * (Math.random() - 0.5);

                final double flowDepthTkh;
                final double tkhUp;
                final double tkhDown;
                switch (kind) {
                case starr:
                    flowDepthTkh = wst - (meanBedHeight + tkh / 100);
                    tkhUp = tkh;
                    tkhDown = 0;
                    break;

                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);

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

    /**
     * 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.
     *
     * Zeitraum Zeitliche Differenz [a]
     * X ≥ 1998 ± 3
     * 1958 ≤ X < 1998 ± 6
     * 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();
        if (soundingYear == null)
            return;

        final int wstYear = waterlevel.getYear();
        if (wstYear < 0)
            return;

        final int maxDifference = getMaxDifferenceYears(soundingYear);

        final int difference = Math.abs(soundingYear - wstYear);
        if (difference > maxDifference) {
            final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.year_difference", null, label, wstYear,
                    soundingYear);
            problems.addProblem(message);
        }
    }

    private int getMaxDifferenceYears(final int year) {

        if (year < 1918)
            return 25;

        if (1918 <= year && year < 1958)
            return 12;

        if (1958 <= year && year < 1998)
            return 6;

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

        final int size = wstKms.size();
        for (int i = 0; i < size - 2; i++) {
            final double kmPrev = wstKms.getKm(i);
            final double kmNext = wstKms.getKm(i + 1);

            /* only check if we are within the calculation range */
            if (calcRange.overlapsRange(new DoubleRange(kmPrev, kmNext))) {
                if (Math.abs(kmPrev - kmNext) > 1) {
                    final String label = wstKms.getName();

                    final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.waterlevel_discretisation", null, label);
                    problems.addProblem(kmPrev, message);
                }
            }
        }
    }

    private BedHeight loadBedHeight(final String soundingId) {

        // REMARK: absolutely unbelievable....
        // The way how bed-heights (and other data too) is accessed is different for nearly every calculation-type
        // throughout flys.
        // The knowledge on how to parse the datacage-ids is spread through the complete code-base...

        // We use here the way on how bed-heights are accessed by the BedDifferenceAccess/BedDifferenceCalculation, but
        // this is plain random
        final String[] parts = soundingId.split(";");

        final BedHeightsArtifact artifact = (BedHeightsArtifact) RiverUtils.getArtifact(parts[0], this.context);

        final Integer bedheightId = artifact.getDataAsInteger("height_id");
        // 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.
        // final String bedheightType = artifact.getDataAsString("type");

        // REMARK: BedDifferences uses this, but we also need the metadata of the BedHeight
        // REMARK: second absolutely awful thing: BedHeight is a hibernate binding class, accessing the database via
        // hibernate stuff
        // 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(double h, double vm, double d50, 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);
    }
}