using System;
using Deltares.Geometry;

namespace Deltares.Stability.Calculation2
{
    public class SpencerCalculation : SlidingCurveCalculation
    {
        public const int CBoundsMin = 0;
        public const int CBoundsInc = 1;
        private const double CDeltaShift = 0.01;
        private const double CInverseSafetyShift = 0.01;
        private const int CMaxIterations = 10;
        private const double CDegreeToRad = Math.PI/180;

        private double delta = 0;
        private double safety;
        private SlidingPlane 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.00001);
            }
        }

        public override double GetSafetyFactor(SlidingCurve slidingCurve)
        {
            if (slidingCurve is SlidingPlane)
            {
                var spencerPlane = (SlidingPlane) slidingCurve;
                this.slidingCurve = spencerPlane;
                // TODO: use sliding circle data instead of internal fields
                return CalculateStabilityFactor(spencerPlane); // to be implemented
            }

            throw new ArgumentException("Sliding curve should be sliding plane");
        }

        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.Spencerplane = fillPlane((SlidingPlane) slidingCurve);

            SlipplaneConstraints constraints = PreProcessCalculation.Constraints;
            sliceCreatorCreator.StriveWidth = slidingCurve.MaximumSliceWidth;
            sliceCreatorCreator.CreateSpencerSlices((SlidingPlane) slidingCurve);
        }

        private void SwapSlices(int firstIndex, int secondIndex, SlidingCurve slidingCurve)
        {
            Slice slice = slidingCurve.Slices[firstIndex];
            slidingCurve.Slices[firstIndex] = slidingCurve.Slices[secondIndex];
            slidingCurve.Slices[secondIndex] = slice;
        }

        private void ReverseData()
        {
            int i;
            int j;
            int Last;
            double Temp;
            var TempPoint = new GeometryPoint();
            TempPoint.X = slidingCurve.LeftPoint.X;
            TempPoint.Z = slidingCurve.LeftPoint.Z;
            slidingCurve.LeftPoint.X = -slidingCurve.RightPoint.X;
            slidingCurve.LeftPoint.Z = -slidingCurve.RightPoint.Z;
            slidingCurve.RightPoint.X = -TempPoint.X;
            slidingCurve.RightPoint.Z = -TempPoint.Z;

            Slice lSlice;
            for (int k = 0; k < slidingCurve.Slices.Count; k++)
            {
                lSlice = slidingCurve.Slices[k];
                // xcoors
                Temp = lSlice.TopLeft.X;
                lSlice.TopLeft.X = -lSlice.TopRight.X;
                lSlice.BottomLeft.X = -lSlice.TopRight.X;
                lSlice.TopRight.X = -Temp;
                lSlice.BottomRight.X = -Temp;

                // ZCoors
                Temp = lSlice.BottomLeft.Y;
                lSlice.BottomLeft.Y = lSlice.BottomRight.Y;
                lSlice.BottomRight.Y = Temp;

                Temp = lSlice.TopLeft.Y;
                lSlice.TopLeft.Y = lSlice.TopRight.Y;
                lSlice.TopRight.Y = Temp;

                Temp = lSlice.LeftForceY;
                lSlice.LeftForceY = lSlice.RightForceY;
                lSlice.RightForceY = Temp;

                Temp = lSlice.LeftForce;
                lSlice.LeftForce = lSlice.RightForce;
                lSlice.RightForce = Temp;

                // Loading
                lSlice.HPoreOnSurface = -lSlice.HPoreOnSurface;
                lSlice.LeftForceAngle = -lSlice.LeftForceAngle;
                // right force angle is missing (is conected to leftforce angle
            }

            // change numbering of slices

            // wissel de slices om (aleen de adressen)
            Last = Convert.ToInt32(slidingCurve.Slices.Count/2);
            for (i = 0; i < Last; i++)
            {
                j = slidingCurve.Slices.Count - 1 - i;
                SwapSlices(i, j, slidingCurve);
            }
        }

        private void DeltaGrens(ref double AMinDelta, ref double AMaxDelta)
        {
            int i;
            double Ltmpphi;
            double LDeltaMin;
            double LDeltaMax;
            AMinDelta = -0.5*Math.PI;
            AMaxDelta = 0.5*Math.PI;
            for (i = 0; i < slidingCurve.Slices.Count; i++)
            {
                Slice LSice = slidingCurve.Slices[i];
                if (LSice.Phi*Constants.CRad < 1.0E-12)
                {
                    Ltmpphi = 1.0E-12;
                }
                else
                {
                    Ltmpphi = LSice.Phi*Constants.CRad;
                }
                // minimale grenswaarde bepalen
                LDeltaMin = -slidingCurve.Slices[i].BottomAngle -
                            (Math.PI + Math.Atan(-safety/Math.Tan(Ltmpphi)));
                AMinDelta = Math.Max(LDeltaMin, AMinDelta);
                // maximale grenswaarde bepalen
                LDeltaMax = LDeltaMin + Math.PI;
                AMaxDelta = Math.Min(LDeltaMax, AMaxDelta);
            }
        }

        private void DetermineResultantMoment()
        {
            int LSliceIndex;
            for (LSliceIndex = 0; LSliceIndex < slidingCurve.Slices.Count; LSliceIndex++)
            {
                slidingCurve.Slices[LSliceIndex].DetermineResultantMoment();
            }
        }

        //=======================================================================================================================
        //Comment: The blue print of the theoretical implications is found in:
        //  https://repos.deltares.nl/repos/pvcsprojects/documents/trunk/DGeoStability/Design/BluePrintGeoStability.doc
        //
        //Date        ID   Modification
        //2010-11-17  Sel  Created
        //=======================================================================================================================}
        private void ForceSolutionWithinRestraints()
        {
            int LSliceIndex;
            double LInverseSafety;
            double LDelta;

            //Confine solution space; the constraints are found in equation (9) 
            LInverseSafety = 1/safety + CInverseSafetyShift;
            if (LInverseSafety < CInverseSafetyShift)
            {
                //  throw EDSpencerPlaneError.Create('New safety factor is negative; no proper solution is anymore expected');
            }

            // The resultant slope range is sought for LInversSafety; then, also valid for LInversSafety - CInverseSafetyShift 
            double LMinimumDelta = double.MinValue;
            double LMaximumDelta = double.MaxValue;

            // Find the Delta range belonging to LSafety 
            for (LSliceIndex = 0; LSliceIndex < slidingCurve.Slices.Count; LSliceIndex++)
            {
                // BottomAngle is defined negative here 
                LMinimumDelta = Math.Max(LMinimumDelta, (-Math.PI/2 - slidingCurve.Slices[LSliceIndex].BottomAngle));
                LMaximumDelta = Math.Min(LMaximumDelta, (Math.PI/2 - Math.Atan(
                    Math.Tan(slidingCurve.Slices[LSliceIndex].Phi*CDegreeToRad)*
                    LInverseSafety) - slidingCurve.Slices[LSliceIndex].BottomAngle));
            }

            // Stay away CDeltaShift from the restraints; reserve space for LDelta - CDeltaShift 
            if (LMinimumDelta > LMaximumDelta - 3*CDeltaShift)
            {
                //       throw EDSpencerPlaneError.Create('No valid range for resultant slope; no proper solution is anymore expected');
            }

            LDelta = Math.Max(LMinimumDelta + 2*CDeltaShift, Math.Min(delta, LMaximumDelta - CDeltaShift));

            //{ Allow adapted proposed shift of resultant slope 
            delta = LDelta;
        }

        private void DetermineInterSliceForces()
            //==============================================
        {
            Slice slice = slidingCurve.Slices[0];
            // The first slice has no left force }
            slice.LeftForce = 0;
            slice.LeftForceY = slice.BottomLeft.Y;
            slice.LeftForceAngle = delta;

            // Determine the right force and slope and pass them to the left values of the next slice }
            for (int sliceIndex = 0; sliceIndex < slidingCurve.Slices.Count; sliceIndex++)
            {
                slice = slidingCurve.Slices[sliceIndex];
                // Determine resultant }
                slice.TmpSafety = safety;
                slice.ResultantAngle = delta;
                slice.DetermineResultantNorm();
                slice.DetermineResultantForce();
            }
            for (int sliceIndex = 0; sliceIndex < slidingCurve.Slices.Count - 1; sliceIndex++)
            {
                slice = slidingCurve.Slices[sliceIndex];
                slice.DetermineInterSliceForces();
                // Pass to next slice }
                slidingCurve.Slices[sliceIndex + 1].LeftForce = slice.RightForce;
                slidingCurve.Slices[sliceIndex + 1].LeftForceY = slice.RightForceY;
                slidingCurve.Slices[sliceIndex + 1].LeftForceAngle = slice.RightForceAngle;
            }

            // It is neat to specify proper end values of the inter-slice slopes, though they are not relevant }
            slice = slidingCurve.Slices[slidingCurve.Slices.Count - 1];
            slice.RightForce = 0;
            slice.RightForceY = slice.BottomRight.Y;
            slice.RightForceAngle = delta;
        }

        private double DetermineForceImbalance()
        {
            //See equation (10) and (13); the right inter-slice force should vanish at the last slice. In (10) the ForceAngle's
            //    are equal. We use this routine for both methods, Constant Resultant Slope and Equally Spaced Action ratio }

            Slice slice = slidingCurve.Slices[slidingCurve.Slices.Count - 1];
            return slice.LeftForce*Math.Cos(slice.LeftForceAngle - slice.RightForceAngle) +
                   slice.ResultantForce*Math.Cos(slice.ResultantAngle - slice.RightForceAngle);
        }

        private double DetermineMomentImbalance()
        {
            // See equation (5); the right inter-slice force should vanish at the last slice }
            Slice slice = slidingCurve.Slices[slidingCurve.Slices.Count - 1];
            return slice.ResultantMoment + slice.LeftForce*Math.Cos(slice.LeftForceAngle)*
                   (slice.LeftForceY - slice.ZCenterBottom - Math.Tan(slice.LeftForceAngle)*slice.Width/2.0);
        }

        private bool IterateConstantResultantSlopes()
        {
            const double CSafetyAccuracy = 0.0001;
            const double CDeltaAccuracy = 0.0001;
            bool result;
            double LRestForce;
            double LRestMoment;
            double LMomentPerSafety;
            double LForcePerSafety;
            double LMomentPerDelta;
            double LForcePerDelta;
            double LInverseSafetyShift;
            double LDeltaShift;

            // Determine present imbalance 
            DetermineInterSliceForces();
            LRestMoment = DetermineMomentImbalance();
            LRestForce = DetermineForceImbalance();

            //  Determine imbalance for variation in resultant slope }
            delta = delta - CDeltaShift;
            DetermineInterSliceForces();
            LMomentPerDelta = (LRestMoment - DetermineMomentImbalance())/CDeltaShift;
            LForcePerDelta = (LRestForce - DetermineForceImbalance())/CDeltaShift;
            delta = delta + CDeltaShift;

            // etermine imbalance for variation in safety factor }
            safety = safety/(1 + safety*CInverseSafetyShift);
            DetermineInterSliceForces();
            LMomentPerSafety = (DetermineMomentImbalance() - LRestMoment)/CInverseSafetyShift;
            LForcePerSafety = (DetermineForceImbalance() - LRestForce)/CInverseSafetyShift;
            safety = safety/(1 - safety*CInverseSafetyShift);

            // ermine the adjustments }
            LDeltaShift = LForcePerSafety*LMomentPerDelta - LForcePerDelta*LMomentPerSafety;
            LInverseSafetyShift = (LMomentPerDelta*LRestForce - LForcePerDelta*LRestMoment)/LDeltaShift;
            LDeltaShift = (LForcePerSafety*LRestMoment - LMomentPerSafety*LRestForce)/LDeltaShift;

            //Check accuracy and adjust }
            result = (Math.Abs(LInverseSafetyShift) < CSafetyAccuracy) && (Math.Abs(LDeltaShift) < CDeltaAccuracy);
            //  Result := ((Abs(LRestForce) < CEpsForce) and (Abs(LRestMoment) < CEpsMom));
            safety = safety/(1 - safety*LInverseSafetyShift);
            delta = delta - LDeltaShift;
            return result;
        }

        private bool IsInterSliceForceOutsideSlice()
        {
            int LSliceIndex;
            for (LSliceIndex = 0; LSliceIndex < slidingCurve.Slices.Count - 1; LSliceIndex++)
            {
                if ((slidingCurve.Slices[LSliceIndex].RightForceY < slidingCurve.Slices[LSliceIndex].BottomRight.Y)
                    || (slidingCurve.Slices[LSliceIndex].RightForceY > slidingCurve.Slices[LSliceIndex].TopRight.Y))
                {
                    return true;
                }
            }
            return false;
        }

        private void GetShearForces()
        {
            Slice LSlice;
            int LSliceIndex;
            double LeftCombination;
            double LRightCombination;
            double LTotalPore;
            double LTotPore;
            double LCosDif;
            double LCosBot;
            double LMass;
            double LNormalStress;
            double LTanPhiBot;
            double LCohBottom;

            for (LSliceIndex = 0; LSliceIndex < slidingCurve.Slices.Count - 1; LSliceIndex++)
            {
                LSlice = slidingCurve.Slices[LSliceIndex];
                if (LSlice.ArcLength > 0)
                {
                    LeftCombination = LSlice.LeftForce*Math.Sin(-LSlice.BottomAngle - LSlice.LeftForceAngle);
                    LRightCombination = LSlice.RightForce*Math.Sin(-LSlice.BottomAngle - LSlice.RightForceAngle);
                    LTotPore =
                        Math.Sqrt(LSlice.VPoreOnSurface*LSlice.VPoreOnSurface +
                                  LSlice.HPoreOnSurface*LSlice.HPoreOnSurface);
                    LTotalPore = LSlice.TotalPorePressure;
                    LCosDif = Math.Cos(LSlice.TopAngle - LSlice.BottomAngle);
                    LCosBot = Math.Cos(-LSlice.BottomAngle);
                    LMass = LSlice.Weight + LSlice.ExternalLoad + LSlice.SigmaSoilQuake*LSlice.Width; // in kN/m
                    LSlice.NormalStress = (LTotPore*LCosDif - LTotalPore*LSlice.ArcLength +
                                           (LMass)*LCosBot - LeftCombination + LRightCombination)/LSlice.ArcLength;
                    LTanPhiBot = Math.Tan(LSlice.Phi*CDegreeToRad);
                    LNormalStress = LSlice.NormalStress;
                    LCohBottom = LSlice.Cohesion;
                    LSlice.ShearStress = (LNormalStress*LTanPhiBot + LCohBottom)/safety;
                }
            }
        }

        /// <summary>
        /// Calculate spencer with
        /// ApplyConstantResultantSlopes 
        /// </summary>
        /// <returns></returns>
        private double CalculateStabilityFactor(SlidingPlane spencerPlane)
        {
            int LIter;
            bool LIsAccurate = false;

            bool isReverseGeometry = slidingCurve.Slices[0].ZCenterTop <
                                     slidingCurve.Slices[slidingCurve.Slices.Count - 1].ZCenterTop;
            // Fill soil,stress data of slices ( in een apparte routine STSliceFill
            if (isReverseGeometry)
            {
                // om uit te kunnene rekenen moet de geometrie gespiegeld worden
                // dit gebeurd door alle xen van een - teken te voorzien en de lamelen
                // andersom te sorteren
                ReverseData();
            }

            //{ Initialize }
            LIter = 0;
            // begin startwaarde stabiliteits factor van 1.1
            safety = Constants.CFMinStart;
            delta = 0;

            DetermineResultantMoment();

            // For safety's sake check restraints }
            ForceSolutionWithinRestraints();
            // Take care of the iteration }
            while (!LIsAccurate && LIter <= CMaxIterations)
            {
                LIsAccurate = IterateConstantResultantSlopes();
                ForceSolutionWithinRestraints();
                LIter++;
                if ((LIter == CMaxIterations) || (safety < 0))
                {
                    /*throw new Exception(
                        "Gradient approach not successful; no proper solution is anymore expected");*/
                    safety = Double.MaxValue;
                    break;
                }
            }

            // Check position of the inter-slice forces }
            if (LIsAccurate)
            {
                DetermineInterSliceForces();
            }

            if (IsInterSliceForceOutsideSlice())
            {
                //excludedSpencerPlanes.Add(spencerPlane);
            }

            if (LIsAccurate)
            {
                GetShearForces();
            }

            if (isReverseGeometry)
            {
                ReverseData();
            }

            spencerPlane.SafetyFactor = safety;

            return safety;
        }

//===========================================================================

        private GeometryPointString fillPlane(SlidingPlane slidingPlane)
        {
            var result = new GeometryPointString();
            for (int i = 0; i < slidingPlane.SpencerSlipPlane.Count; i++)
            {
                var LPoint = new GeometryPoint();
                LPoint.X = slidingPlane.SpencerSlipPlane[i].X;
                LPoint.Z = slidingPlane.SpencerSlipPlane[i].Y;
                result.Points.Add(LPoint);
            }
            return result;
        }
    }
}