using System;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Soils;
using Deltares.Standard;
using Deltares.Standard.Logging;

namespace Deltares.Stability.Calculation2
{
    public class BishopCalculation : SlidingCurveCalculation
    {
        public ModelOptions ModelOption { get; set; }

        public override void CreateSlices(SlidingCurve slidingCurve)
        {
            var sliceCreatorCreator = new SlicesCreator();
            sliceCreatorCreator.ChangingPoints = PreProcessCalculation.ChangingPoints;
            sliceCreatorCreator.SetSoilProfile2D(PreProcessCalculation.Profile2D);
            sliceCreatorCreator.SetWaternet(PreProcessCalculation.Waternet);
            sliceCreatorCreator.SetSurfaceLine(PreProcessCalculation.SurfaceLine);

            sliceCreatorCreator.StriveWidth = slidingCurve.MaximumSliceWidth;
            sliceCreatorCreator.CreateBishopSlices((SlidingCircle) slidingCurve);
        }

        public override void FillSlices(SlidingCurve slidingCurve)
        {
            var fillSlicesFiller = new SlicesFiller();
            fillSlicesFiller.PreProces = PreProcessCalculation;
            for (int i = 0; i < slidingCurve.Slices.Count; i++)
            {
                fillSlicesFiller.FillSlicesByInterpolation(slidingCurve.Slices[i], 0.0);
            }
        }

        public override double GetSafetyFactor(SlidingCurve slidingCurve)
        {
            if (slidingCurve is SlidingCircle)
            {
                var slidingCircle = (SlidingCircle) slidingCurve;
                return CalculateStabilityFactor(slidingCircle); // to be implemented
            }
            else
            {
                throw new ArgumentException("Sliding curve should be circle");
            }
        }

        /// <summary>
        /// Determines the soil moment due to weight of the soil
        /// </summary>
        /// <returns></returns>
        private double DetermineDrivingSoilMoment(SlidingCircle slidingCircle)
        {
            double moment = 0;
            for (int i = 0; i < slidingCircle.Slices.Count; i++)
            {
                double HorDistance = slidingCircle.Slices[i].XCenter - slidingCircle.Center.X;
                // Bereken aandrijvend moment : dx * lamel gewicht
                moment = moment - HorDistance*slidingCircle.Slices[i].Weight;
            }

            return moment;
        }

        private double DetermineWaterMoment(SlidingCircle slidingCircle)
        {
            double moment = 0;
            for (int i = 0; i < slidingCircle.Slices.Count; i++)
            {
                // Calculate watermoment due to waterforces on slidingCurve.Slices
                if ((Math.Abs(slidingCircle.Slices[i].VPoreOnSurface - 0.0) > 0.001) || (Math.Abs(slidingCircle.Slices[i].HPoreOnSurface - 0.0) > 0.001))
                {
                    double horDistance = slidingCircle.Slices[i].XCenter - slidingCircle.Center.X;
                    double vertDistance = slidingCircle.Center.Y - slidingCircle.Slices[i].ZCenterTop;
                    moment = moment + horDistance*slidingCircle.Slices[i].VPoreOnSurface + vertDistance*slidingCircle.Slices[i].HPoreOnSurface;
                }
            }
            var circle = new Circle()
                {
                    XleftSide = slidingCircle.LeftPoint.X,
                    XrightSide = slidingCircle.RightPoint.X,
                    Xcentre = slidingCircle.Center.X,
                    Zcentre = slidingCircle.Center.Y,
                    Radius = slidingCircle.Radius
                };
            moment = moment + WaterMomentOnVerticals(circle);
            return moment;
        }

        private double BishopMethode(SlidingCircle slidingCircle)
        {
            double result;
            const double Pir = Math.PI/180.0;
            int i;
            double Fn;
            double Fo;
            bool Ready;
            double Salfa;
            double TgPhi;
//          double LSinPhi;
//          double LCosPhi;
//          double LCosDilet;
//          double LSinDilet;
            double LTau;
            double Coh;
            double alfa;
            double h;
            double tgalfa;
            //double LCohTerm;
            //double LSigTerm;

            int numberOfIterations = 0;

            // Initialisatie
            slidingCircle.ResistingMoment = 0.0;
            Fn = StartValueSafetyfactor;
            Ready = false;

            slidingCircle.DrivingMoment = slidingCircle.WaterMoment + slidingCircle.SoilMoment + slidingCircle.LoadMoment;

            // Controleer aandrijvend moment op 1 kNm
            if (Math.Abs(slidingCircle.DrivingMoment) < 0.001)
            {
                Ready = true;
            }

            // Start iteratie
            if (Ready == false)
            {
                do
                {
                    numberOfIterations++;
                    Fo = Fn;
                    slidingCircle.ResistingMoment = 0.0;
                    for (i = 0; i < slidingCircle.Slices.Count; i++)
                    {
                        Slice slice = slidingCircle.Slices[i];
                        Coh = slice.Cohesion;
                        // cohesie
//                        LSinPhi = Math.Sin(slice.Phi * Pir);
//                        LCosPhi = System.Math.Cos(slice.Phi * Pir);
                        //                    LCosDilet = Math.Cos(slice.Dilatancy * Pir);
                        //LSinDilet = Math.Sin(slice.Dilatancy * Pir);
                        TgPhi = Math.Tan(slice.Phi*Pir);
                        // formula based on Tau-alfa =
                        // 1/Fo * (coh * cotgphi + Salfa')* cospsi*sinphi/(1 - sinphi*sinpsi)
                        // (after tuenissen))
                        // note that if psi = phi  then  tau = coh + salfa * tgphi
                        //                       LSigTerm = LCosDilet * LSinPhi / (1 - LSinDilet * LSinPhi);
//                        LCohTerm = LCosDilet * LCosPhi / (1 - LSinDilet * LSinPhi);
                        // Tau-alfa = Coh*LCohTerm + Salfa'*LSigTerm
                        alfa = slice.BottomAngle;
                        // Hoek onderkant lamel [rad]
                        if (slidingCircle.DrivingMoment > 0.0)
                        {
                            alfa = -alfa;
                            // Passieve alfa nagatief maken
                        }
                        // Afsnuiten indien nodig
                        h = 0.5*Math.Atan(TgPhi/Fo) - 0.25*Math.PI;
                        if (alfa <= h)
                        {
                            alfa = h;
                        }
                        tgalfa = Math.Tan(alfa);

                        // Bereken het tegen werkend  moment
                        //                       LTau = Coh * LCohTerm + slice.EffectiveStress * LSigTerm;
                        LTau = Coh + slice.EffectiveStress*TgPhi;
                        slidingCircle.ResistingMoment = slidingCircle.ResistingMoment + slice.ArcLength*LTau/(1 + tgalfa*TgPhi/Fo);

                        // Bereken het tegen werkend  moment
                        // FResistingMom := FResistingMom + ArcLength * (coh + EffectiveStress * tgphi) /
                        // (1 + tgalfa * tgphi / Fo);
                        // Vertikaal Evenwicht   : Sigma-v' = Salfa' + Tau-alfa*Tgalfa
                        // en Tau volgens Bishop : Tau-alfa = 1/Fo * (Coh*LCohTerm + Salfa'*LSigTerm
                        // Uit deze twee vergelijkingen kan Salfa opgelost worden  :
                        Salfa = (Fo*slice.EffectiveStress - Coh*tgalfa)/(Fo + TgPhi*tgalfa);
                        // Salfa := (Fo * EffectiveStress - LLocCoh * tgalfa) / (Fo + tgalfa * tgphi);
                        // Nieuwe Tau - alfa
                        slice.NormalStress = Salfa;
                        slice.ShearStress = (Coh + Salfa*TgPhi)/Fo;
                    }

                    slidingCircle.ResistingMoment = slidingCircle.ResistingMoment*slidingCircle.Radius;
                    Fn = slidingCircle.ResistingMoment/Math.Abs(slidingCircle.DrivingMoment);
                    // Sla ongeitereerd moment, dus bij 1e iteratie op, om te kunnen uitvoeren
                    if (numberOfIterations == 1)
                    {
                        slidingCircle.ResistingMomentUniterated = slidingCircle.ResistingMoment;
                    }
                    // Stoppen als veiligheidscoefficient oud en nieuw ongeveer gelijk zijn
                    if (Math.Abs(Fn - Fo) < 0.001)
                    {
                        Ready = true;
                    }

                    // Teveel iteratie dan stoppen
                    if (numberOfIterations > Constants.CMaxIter)
                    {
                        Ready = true;
                    }
                    // Als tegenwerkend moment = 0, dan ook stoppen.
                    if ((Math.Abs(slidingCircle.ResistingMoment) < 0.001) || (Fn < 0.001))
                    {
                        Ready = true;
                    }
                } while (!(Ready));
            }

            result = Fn;
            if (Math.Abs(slidingCircle.DrivingMoment) < 0.001)
            {
                result = 1001;
                LogManager.Add(new LogMessage(LogMessageType.Warning, this, "Driving moment too small"));
            }
            else if (numberOfIterations >= Constants.CMaxIter)
            {
                result = 1001;
                LogManager.Add(new LogMessage(LogMessageType.Warning, this, "More than max " + Constants.CMaxIter + " iterations"));
            }
            else if (Math.Abs(slidingCircle.ResistingMoment) < 0.001 || Fn < 0.001)
            {
                result = 1001;
                LogManager.Add(new LogMessage(LogMessageType.Warning, this, "Resisting moment too small"));
            }

            return result;
        }

        /// <summary>
        /// Calculates the stability calculation for Bishop method
        /// </summary>
        /// <returns></returns>
        private double CalculateStabilityFactor(SlidingCircle slidingCircle)
        {
            slidingCircle.SoilMoment = DetermineDrivingSoilMoment(slidingCircle);
            if (PreProcessCalculation.Waternet.PhreaticLine != null)
            {
                slidingCircle.WaterMoment = DetermineWaterMoment(slidingCircle);
            }

            slidingCircle.LoadMoment = 0;
            foreach (var uniformLoad in PreProcessCalculation.UniformLoads)
            {
                slidingCircle.LoadMoment += uniformLoad.GetMomentRoundX(slidingCircle.Center.X,
                                                                        slidingCircle.LeftPoint.X,
                                                                        slidingCircle.RightPoint.X);
            }

            //
            // Nu aandrijvende momenten en lamel grootheden bekend zijn,
            // kunnen de tegenwerkende momenten en de stabiliteits factor bepaald worden
            //

            slidingCircle.SafetyFactor = BishopMethode(slidingCircle);

            return slidingCircle.SafetyFactor;
        }
    }
}