using System;
using Deltares.Mathematics;

namespace Deltares.Dam.Data
{
    public class DesignPointCalculation
    {
        public enum ExceedingSet
        {
            tenThousend,
            fourThousend,
            twoThousend,
            TwelfFifty
        };

        private double FAlfaH;
        private double FBeta;
        private double[] FBetas;
        private double FDecimate;
        private double FDelta, FDeltaFactor;
        private string FErrorMessage;
        private ExceedingSet FExceed;
        private double FHOntwerp, FHwaterAlfa;
        private bool FIsMaxLevelUsed;
        private int FLevelCount;
        private double[] FLocBetas, FLocLevels;
        private int FLoopCount;
        private double FMHW;
        private double FMaxLevel;
        private double[] FRepetitiveBeta;
        private int FRepetitiveCounter;
        private double[] FRepetitiveHNew;
        private double[] FWaterlevels;

        public double Beta
        {
            get
            {
                return FBeta;
            }
            set
            {
                FBeta = value;
            }
        }

        public double HOntwerp
        {
            get
            {
                return FHOntwerp;
            }
            set
            {
                FHOntwerp = value;
            }
        }

        public double HwaterAlfa
        {
            get
            {
                return FHwaterAlfa;
            }
            set
            {
                FHwaterAlfa = value;
            }
        }

        public double AlfaH
        {
            get
            {
                return FAlfaH;
            }
            set
            {
                FAlfaH = value;
            }
        }

        public double MHW
        {
            get
            {
                return FMHW;
            }
            set
            {
                FMHW = value;
            }
        }

        public double Decimate
        {
            get
            {
                return FDecimate;
            }
            set
            {
                FDecimate = value;
            }
        }

        public double MaxLevel
        {
            get
            {
                return FMaxLevel;
            }
            set
            {
                FMaxLevel = value;
            }
        }

        public int LoopCount
        {
            get
            {
                return FLoopCount;
            }
            set
            {
                FLoopCount = value;
            }
        }

        public string ErrorMessage
        {
            get
            {
                return FErrorMessage;
            }
            set
            {
                FErrorMessage = value;
            }
        }

        public ExceedingSet Exceed
        {
            get
            {
                return FExceed;
            }
            set
            {
                FExceed = value;
            }
        }

        public bool IsMaxLevelUsed
        {
            get
            {
                return FIsMaxLevelUsed;
            }
            set
            {
                FIsMaxLevelUsed = value;
            }
        }

        public double[] Waterlevels
        {
            get
            {
                return FWaterlevels;
            }
            set
            {
                FWaterlevels = value;
            }
        }

        public double[] Betas
        {
            get
            {
                return FBetas;
            }
            set
            {
                FBetas = value;
            }
        }

        public bool CalculateTheWaterDesignpoint()
        {
            bool LReady;
            double UGumbel;
            double dBetadH;
            double HNew;
            double Gamma;
            double Inverse;
            double Ehn;
            double Ez;
            double Shn;
            double Sz;
            double axu;
            double alfa;
            double CMaxDouble = 1.7E+308;

            LReady = false;
            FDeltaFactor = 1;
            HNew = 0;

            // determine the parameters of the equivalent normal distribution
            // fill and sort the local beta array
            FillLocalBetas();
            SortLocalBetas();
            FRepetitiveCounter = 0;
            FRepetitiveBeta = new double[30];
            FRepetitiveHNew = new double[30];

            if (FIsMaxLevelUsed)
            {
                AddMaxlevelToLocalDeltas();
            }
            if (FDecimate > 0.001)
            {
                UGumbel = GetGumbel();
                alfa = 2.3/FDecimate;
                // Startwaarde voor H
                FHOntwerp = FMHW;
                FLoopCount = 0;

                do
                {
                    // find beta on curve
                    FBeta = GetBeta(FHOntwerp);
                    axu = -alfa*(FHOntwerp - UGumbel);
                    Inverse = GetInverse(axu);
                    if (Math.Abs(Inverse) < (Math.Sqrt(CMaxDouble) - 1))
                    {
                        Gamma = GetGamma(alfa, axu, Inverse);
                        Ehn = HOntwerp - Inverse*Gamma*FDecimate;
                        Shn = Gamma*FDecimate;
                        dBetadH = GetSlope(HOntwerp);
                        FLoopCount++;
                        if (FLoopCount >= 20)
                        {
                            LReady = IsNonIterative(FBeta, HNew);
                        }
                        if (!LReady)
                        {
                            Ez = FBeta + dBetadH*(Ehn - HOntwerp);
                            Sz = Math.Sqrt(1 + dBetadH*dBetadH*Shn*Shn);
                            FBeta = Ez/Sz;
                            FAlfaH = dBetadH*Shn/Sz;
                            HNew = UGumbel -
                                   FDecimate/2.3*
                                   Math.Log(-Math.Log(Probabilistic.Probabilistic.NormalDistribution(-FBeta*FAlfaH), Math.E), Math.E);
                            FDelta = HNew - FHOntwerp;
                            LReady = (Math.Abs(FDelta) < 0.001);
                            FHOntwerp = FHOntwerp + FDeltaFactor*FDelta;
                        }
                    }
                    else
                    {
                        //error find old FHOntwerp and use smaller steps 
                        //only if 1 step is already calculated
                        FDeltaFactor = FDeltaFactor*0.5;
                        if (FLoopCount > 0)
                        {
                            FHOntwerp = HNew - FDelta;
                            FDelta = 0.5*FDelta;
                            FHOntwerp = FHOntwerp + FDeltaFactor*FDelta;
                            FLoopCount++;
                        }
                        else
                        {
                            var designPointCalculationException = new Exception("Eternal loop");
                            throw designPointCalculationException;
                        }
                    }
                } while (!LReady && (FLoopCount <= 50));

                if (LReady)
                {
                    FHwaterAlfa = GetNearestLevel(FHOntwerp);
                }
                BetaOutputWaterlevelsIncluded(LReady);
            }
            return LReady;
        }

        private void FillLocalBetas()
        {
            int i;

            FLevelCount = FBetas.GetLength(0);
            FLocBetas = new double[FLevelCount + 1]; // is cheating a bit, should be a dynamic array to accomodate AddMaxlevelToLocalDeltas
            FLocLevels = new double[FLevelCount + 1]; // is cheating a bit, should be a dynamic array to accomodate AddMaxlevelToLocalDeltas
            for (i = FBetas.GetLowerBound(0); i <= FBetas.GetUpperBound(0); i++)
            {
                // put waterlevels and betas in local arrays
                FLocLevels[i] = FWaterlevels[i];
                FLocBetas[i] = FBetas[i];
            }
        }

        private void SwapLevels(int FromLevel, int ToLevel)
        {
            double Ltmp1;
            double Ltmp2;

            Ltmp1 = FLocLevels[ToLevel];
            Ltmp2 = FLocBetas[ToLevel];
            FLocLevels[ToLevel] = FLocLevels[FromLevel];
            FLocBetas[ToLevel] = FLocBetas[FromLevel];
            FLocLevels[FromLevel] = Ltmp1;
            FLocBetas[FromLevel] = Ltmp2;
        }

        private void SortLocalBetas()
        {
            int i;
            int j;

            for (i = 0; i <= FLevelCount - 1; i++)
            {
                for (j = i; j <= FLevelCount - 1; j++)
                {
                    if (FLocLevels[j] < FLocLevels[i])
                    {
                        SwapLevels(j, i);
                    }
                }
            }
        }

        private void AddMaxlevelToLocalDeltas()
        {
            int i;
            int LNum;
            double LX1;
            double LY1;
            double LMaxBeta;
            double[] LLevel;
            double[] LBeta;

            // if maxlevel lower than lowest level present extrapolate lower end create horizontal part 
            if (FMaxLevel <= FLocLevels[0])
            {
                LMaxBeta = Routines2D.LinInpolY(FLocLevels[0], FLocBetas[0], FLocLevels[1], FLocBetas[1], FMaxLevel);
                //setlength(LLevel, 3);
                LLevel = new double[3];
                //setlength(LBeta, 3);
                LBeta = new double[3];
                LX1 = FMaxLevel - 1;
                LY1 = Routines2D.LinInpolY(FLocLevels[0], FLocBetas[0], FLocLevels[1], FLocBetas[1], LX1);

                LLevel[0] = LX1;
                LBeta[0] = LY1;

                LLevel[1] = FMaxLevel;
                LBeta[1] = LMaxBeta;

                LLevel[2] = FMaxLevel + 1;
                LBeta[2] = LMaxBeta;

                FLevelCount = 3;
                //Setlength(FLocBetas, FLevelCount);
                //Setlength(FLocLevels, FLevelCount);
                for (i = 0; i <= 2; i++)
                {
                    FLocBetas[i] = LBeta[i];
                    FLocLevels[i] = LLevel[i];
                }
            }
            else
            {
                if (FMaxLevel > FLocLevels[FLevelCount - 1])
                {
                    //if maxlevel higher than highestt level present then add a horizontal part after maxlevel
                    LMaxBeta = Routines2D.LinInpolY(FLocLevels[FLevelCount - 1], FLocBetas[FLevelCount - 1],
                                                    FLocLevels[FLevelCount - 2], FLocBetas[FLevelCount - 2], FMaxLevel);
                    // Te highest point is in levelcount - 1 
                    FLocLevels[FLevelCount - 1] = FMaxLevel;
                    FLocBetas[FLevelCount - 1] = LMaxBeta;
                    // Add a horizontal part
                    FLevelCount++;
                    //Setlength(FLocBetas, FLevelCount);
                    //Setlength(FLocLevels, FLevelCount);
                    FLocLevels[FLevelCount - 1] = FMaxLevel + 1;
                    FLocBetas[FLevelCount - 1] = LMaxBeta;
                }
                else
                {
                    LNum = 0;
                    LMaxBeta = 0;
                    // GeometryPoint is somewhere in between
                    for (i = 0; i <= FLevelCount - 2; i++)
                    {
                        if ((FMaxLevel >= FLocLevels[i]) && (FMaxLevel <= FLocLevels[i + 1]))
                        {
                            LMaxBeta = Routines2D.LinInpolY(FLocLevels[i], FLocBetas[i], FLocLevels[i + 1],
                                                            FLocBetas[i + 1], FMaxLevel);
                            LNum = i + 1;
                            FLocLevels[LNum] = FMaxLevel;
                            FLocBetas[LNum] = LMaxBeta;
                            break;
                        }
                    }

                    if (LNum < FLevelCount - 1)
                    {
                        for (i = LNum; i <= FLevelCount - 1; i++)
                        {
                            FLocBetas[i] = LMaxBeta;
                        }
                    }
                    else
                    {
                        // LNum is levelcount - 1
                        // Add horizontal part
                        FLevelCount++;
                        FLocLevels[FLevelCount - 1] = FMaxLevel + 1;
                        FLocBetas[FLevelCount - 1] = LMaxBeta;
                    }
                }
            }
        }

        private double GetGumbel()
        {
            double LFaal;
            double result;

            LFaal = 1.0/1250;
            switch (FExceed)
            {
                case ExceedingSet.tenThousend:
                    LFaal = 1.0/10000;
                    break;
                case ExceedingSet.fourThousend:
                    LFaal = 1.0/4000;
                    break;
                case ExceedingSet.twoThousend:
                    LFaal = 1.0/2000;
                    break;
            }
            result = Math.Log(-Math.Log(1 - LFaal, Math.E), Math.E)/(2.3/FDecimate) + FMHW;
            return result;
        }

        private double GetBeta(double HOntwerp)
        {
            int i;
            bool Ready;
            double Result;
            double CEpsMin = 0.0001;

            i = -1;
            Result = 0;

            //... Is HOntwerp on a point ...
            do
            {
                i++;
                Ready = (Math.Abs(HOntwerp - FLocLevels[i]) < CEpsMin);
                if (Ready)
                {
                    Result = FLocBetas[i];
                }
            } while ((!Ready) && ((i < FLevelCount - 1)));

            if (!Ready)
            {
                //... Is HOntwerp between the points ...
                i = -1;
                do
                {
                    i++;
                    Ready = (HOntwerp < Math.Max(FLocLevels[i], FLocLevels[i + 1])) &&
                            (HOntwerp > Math.Min(FLocLevels[i], FLocLevels[i + 1]));
                    if (Ready)
                    {
                        // beta in between
                        Result = Routines2D.LinInpolY(FLocLevels[i], FLocBetas[i], FLocLevels[i + 1], FLocBetas[i + 1],
                                                      HOntwerp);
                    }
                } while ((!Ready) && (i < FLevelCount - 2));

                if (!Ready)
                {
                    //... HOntwerp is outside the points range...
                    if (HOntwerp > FLocLevels[FLevelCount - 1])
                    {
                        Result = Routines2D.LinInpolY(FLocLevels[FLevelCount - 2], FLocBetas[FLevelCount - 2],
                                                      FLocLevels[FLevelCount - 1],
                                                      FLocBetas[FLevelCount - 1], HOntwerp);
                    }
                    else
                    {
                        if (FLevelCount > 1)
                        {
                            Result = Routines2D.LinInpolY(FLocLevels[0], FLocBetas[0], FLocLevels[1], FLocBetas[1],
                                                          HOntwerp);
                        }
                        else
                        {
                            Result = FLocBetas[0];
                        }
                    }
                }
            }
            return Result;
        }

        private double GetInverse(double axu)
        {
            double value;
            double lresult;

            value = Math.Exp(-Math.Exp(axu));
            lresult = Probabilistic.Probabilistic.NormalDistrInverse(value);
            return lresult;
        }

        private double GetGamma(double alfa, double axu, double Inverse)
        {
            double pdf;
            double LResult;
            pdf = alfa*Math.Exp(axu - Math.Exp(axu));
            LResult = Math.Exp(-0.5*Inverse*Inverse)/(pdf*Math.Sqrt(2*Math.PI)*FDecimate);
            return LResult;
        }

        private double Slope(double L1, double B1, double L2, double B2)
        {
            double cTiny = 1.0e-8;

            // If the levels (L) are the same the beta's should be the same as well
            if ((Math.Abs(B2 - B1) < cTiny) || (Math.Abs(L2 - L1) < cTiny))
            {
                return 0;
            }
            else
            {
                return (B2 - B1)/(L2 - L1);
            }
        }

        private double GetSlope(double heigth)
        {
            int i;
            double LResult;

            LResult = 0;
            if (heigth <= FLocLevels[1])
            {
                LResult = Slope(FLocLevels[0], FLocBetas[0], FLocLevels[1], FLocBetas[1]);
            }
            else
            {
                if (heigth >= FLocLevels[FLevelCount - 2])
                {
                    LResult = Slope(FLocLevels[FLevelCount - 2], FLocBetas[FLevelCount - 2],
                                    FLocLevels[FLevelCount - 1], FLocBetas[FLevelCount - 1]);
                }
                else
                {
                    for (i = 1; i <= FLevelCount - 3; i++)
                    {
                        if ((heigth >= FLocLevels[i]) && (heigth <= FLocLevels[i + 1]))
                        {
                            LResult = Slope(FLocLevels[i], FLocBetas[i], FLocLevels[i + 1], FLocBetas[i + 1]);
                        }
                    }
                }
            }
            return LResult;
        }

        private double GetNearestLevel(double HOntwerp)
        {
            int i;
            double value;
            bool ready;

            ready = false;
            value = FLocLevels[0];
            if (HOntwerp <= FLocLevels[0])
            {
                value = FLocLevels[0];
            }
            else
            {
                if (HOntwerp >= FLocLevels[FLevelCount - 1])
                {
                    value = FLocLevels[FLevelCount - 1];
                }
                else
                {
                    i = -1;
                    do
                    {
                        i++;
                        if ((HOntwerp > FLocLevels[i]) && (HOntwerp < FLocLevels[i + 1]))
                        {
                            ready = true;
                            if (Math.Abs(HOntwerp - FLocLevels[i]) < Math.Abs(HOntwerp - FLocLevels[i + 1]))
                            {
                                value = FLocLevels[i];
                            }
                            else
                            {
                                value = FLocLevels[i + 1];
                            }
                        }
                    } while (!ready && (i < FLevelCount - 2));
                }
            }
            return value;
        }

        private void BetaOutputWaterlevelsIncluded(bool Ready)
        {
            FErrorMessage = "";
            if (!Ready)
            {
                FErrorMessage = " No convergence after " + FLoopCount.ToString() + " steps in the iteration process";
            }
        }

        private bool IsRepetitiveHNew(double hNew)
        {
            int repCount = 0;
            for (int i = 0; i < FRepetitiveCounter; i++)
            {
                if (Math.Abs(FRepetitiveHNew[i] - hNew) < 0.0001)
                {
                    repCount ++;
                }
            }
            return (repCount >= 4);
        }

        private double GetAverageBeta()
        {
            double average = 0;
            for (int i = 0; i < FRepetitiveCounter; i++)
            {
                average = average + FRepetitiveBeta[i];
            }
            return average/(FRepetitiveCounter);
        }

        private bool IsNonIterative(double beta, double hNew)
        {
            if (IsRepetitiveHNew(hNew))
            {
                FBeta = GetAverageBeta();
                return true;
            }
            else
            {
                FRepetitiveBeta[FRepetitiveCounter] = beta;
                FRepetitiveHNew[FRepetitiveCounter] = hNew;
                FRepetitiveCounter++;
                return false;
            }
        }
    }
}