using System;
using Deltares.Standard;
using Deltares.Standard.Logging;

namespace Deltares.Stability.Calculation2
{
    public class UpliftVanCalculation : SlidingCurveCalculation
    {
        // internal data
        private double activeDrivingDMoment = 0;
        private double activeDrivingMoment = 0;
        private double activeForce = 0;
        private double activeForce0 = 0;
        private double activeResistingMoment = 0;
        private double drivingMoment = 0;
        private double horForce = 0;
        private double horForce0 = 0;
        private double loadMomentLeft = 0;
        private double loadMomentRight = 0;
        private double passiveDrivingMoment = 0;

        private double passiveForce = 0;
        private double passiveForce0 = 0;
        private double passiveResistingMoment = 0;
        private double radLeft = 0;
        private double radRight = 0;
        private double resistingMom = 0;
        private double resistingMom0 = 0;
        private bool reverseGeom = false;
        private double soilMomentLeft = 0;
        private double soilMomentRight = 0;
        private double totDrivingLeft = 0;
        private double totDrivingRight = 0;
        private double watermomentLeft = 0;
        private double watermomentRight = 0;
        private double xActiveCentre = 0;

        // centre point data in geometry
        private double xMidLeft = 0;
        private double xMidRight = 0;
        private double xPassiveCentre = 0;
        private double zActiveCentre = 0;
        private double zMidLeft = 0;
        private double zMidRight = 0;
        private double zPassiveCentre = 0;
        private double zTangent = 0;

        // destructor Destroy; override;
        // properties
        public double XActiveCentre
        {
            get
            {
                return xActiveCentre;
            }
            set
            {
                xActiveCentre = value;
            }
        }

        public double ZActiveCentre
        {
            get
            {
                return zActiveCentre;
            }
            set
            {
                zActiveCentre = value;
            }
        }

        public double XPassiveCentre
        {
            get
            {
                return xPassiveCentre;
            }
            set
            {
                xPassiveCentre = value;
            }
        }

        public double ZPassiveCentre
        {
            get
            {
                return zPassiveCentre;
            }
            set
            {
                zPassiveCentre = value;
            }
        }

        public double ActiveRadius
        {
            get
            {
                return GetActiveRadius();
            }
        }

        public double PassiveRadius
        {
            get
            {
                return GetPassiveRadius();
            }
        }

        public double ZTangent
        {
            get
            {
                return zTangent;
            }
            set
            {
                zTangent = value;
            }
        }

        public double ActiveDrivingMoment
        {
            get
            {
                return activeDrivingMoment;
            }
            set
            {
                activeDrivingMoment = value;
            }
        }

        public double PassiveDrivingMoment
        {
            get
            {
                return passiveDrivingMoment;
            }
            set
            {
                passiveDrivingMoment = value;
            }
        }

        public double ActiveResistingMoment
        {
            get
            {
                return activeResistingMoment;
            }
            set
            {
                activeResistingMoment = value;
            }
        }

        public double PassiveResistingMoment
        {
            get
            {
                return passiveResistingMoment;
            }
            set
            {
                passiveResistingMoment = value;
            }
        }

        public double ResistingMom0
        {
            get
            {
                return resistingMom0;
            }
            set
            {
                resistingMom0 = value;
            }
        }

        public double ResistingMom
        {
            get
            {
                return resistingMom;
            }
            set
            {
                resistingMom = value;
            }
        }

        public double DrivingMoment
        {
            get
            {
                return drivingMoment;
            }
            set
            {
                drivingMoment = value;
            }
        }

        public double ActiveForce
        {
            get
            {
                return activeForce;
            }
        }

        public double PassiveForce
        {
            get
            {
                return passiveForce;
            }
        }

        public double HorForce
        {
            get
            {
                return horForce;
            }
        }

        public double ActiveForce0
        {
            get
            {
                return activeForce0;
            }
        }

        public double PassiveForce0
        {
            get
            {
                return passiveForce0;
            }
        }

        public double HorForce0
        {
            get
            {
                return horForce0;
            }
        }

        public double WatermomentLeft
        {
            get
            {
                return watermomentLeft;
            }
            set
            {
                watermomentLeft = value;
            }
        }

        public double WatermomentRight
        {
            get
            {
                return watermomentRight;
            }
            set
            {
                watermomentRight = value;
            }
        }

        public double LoadMomentRight
        {
            get
            {
                return loadMomentRight;
            }
            set
            {
                loadMomentRight = value;
            }
        }

        public double LoadMomentLeft
        {
            get
            {
                return loadMomentLeft;
            }
            set
            {
                loadMomentLeft = value;
            }
        }

        public double DetermineWaterMomentOnVertical(double poreTop, double poreBottom, double zTop, double zBottom, double zCentre)
        {
            double lDZ = Math.Abs(zTop - zBottom);
            double lRectPore = Math.Min(poreTop, poreBottom);
            double lForce = lRectPore*lDZ;
            double lWaterMoment;
            lWaterMoment = lForce*(zCentre - 0.5*(zTop + zBottom));
            // triangulair part if present
            lForce = 0.5*lDZ*Math.Abs(poreTop - poreBottom);

            double lArm;
            if ((poreTop < poreBottom))
            {
                lArm = (zCentre - zBottom) - lDZ/3;
            }
            else
            {
                lArm = (zCentre - zTop) + lDZ/3;
            }
            lWaterMoment = lWaterMoment + lForce*lArm;
            return lWaterMoment;
        }

        public override void FillSlices(SlidingCurve slidingCurve)
        {
            SlidingDualCircle slidingDualCurve = (SlidingDualCircle) slidingCurve;

            PrepareCalculation(slidingDualCurve);

            double deltaValue = 0;

            for (int i = 0; i < slidingCurve.Slices.Count; i++)
            {
                var sliceFillerFiller = new SlicesFiller();
                sliceFillerFiller.PreProces = PreProcessCalculation;
                if (i < slidingDualCurve.RightNumber && i > slidingDualCurve.LeftNumber)
                {
                    // in case of Uplift calculation the most unfavourable plane is searched for,
                    // a plane of 2 cm  above the slide plane is compared with a plane 2 cm below the slide plane
                    deltaValue = 0.02;
                    sliceFillerFiller.FillSlicesByInterpolation(slidingCurve.Slices[i], deltaValue);
                    var slice = new Slice();

                    deltaValue = -0.02;
                    sliceFillerFiller.FillSlicesByInterpolation(slice, deltaValue);

                    // Todo AANSLUITEN
                    /*
                    if ((DetermineShear(Slices[i]) > DetermineShear(LSlice)))
                    {
                        Slices[i].CopySliceData(LSlice);
                    }
                     */
                }
                else
                {
                    deltaValue = 0;
                    sliceFillerFiller.FillSlicesByInterpolation(slidingCurve.Slices[i], deltaValue);
                }
            }
        }

        public override double GetSafetyFactor(SlidingCurve slidingCurve)
        {
            if (slidingCurve is SlidingDualCircle)
            {
                var slidingCircle = (SlidingDualCircle) slidingCurve;

                // TODO: use sliding circle data instead of internal fields
                return GetStabilityFactor(slidingCircle); // to be implemented
            }
            else
            {
                throw new ArgumentException("Sliding curve should be dual circle");
            }
        }

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

        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.CreateUpliftVanSlices((SlidingDualCircle) slidingCurve);
        }

        /// <summary>
        /// 
        /// Determine left and right centrepoints
        /// This is more explicit than active and 
        /// passive centrepoints
        /// 
        /// </summary>
        private void PrepareCalculation(SlidingDualCircle slidingCurve)
        {
            xPassiveCentre = slidingCurve.PassiveCircle.X;
            zPassiveCentre = slidingCurve.PassiveCircle.Y;
            xActiveCentre = slidingCurve.ActiveCircle.X;
            zActiveCentre = slidingCurve.ActiveCircle.Y;
            zTangent = slidingCurve.TangentLine;

            if ((xActiveCentre < xPassiveCentre + Constants.CAlmostZero))
            {
                xMidLeft = xActiveCentre;
                zMidLeft = zActiveCentre;
                radLeft = zActiveCentre - zTangent;
                xMidRight = xPassiveCentre;
                zMidRight = zPassiveCentre;
                radRight = zPassiveCentre - zTangent;
            }
            else
            {
                xMidLeft = xPassiveCentre;
                zMidLeft = zPassiveCentre;
                radLeft = zPassiveCentre - zTangent;
                xMidRight = xActiveCentre;
                zMidRight = zActiveCentre;
                radRight = zActiveCentre - zTangent;
            }
        }

        private double DrivingSoilMomentLeft(SlidingDualCircle slidingCurve)
        {
            double moment = 0;

            for (int i = 0; i <= slidingCurve.LeftNumber; i++)
            {
                double horDistance = xMidLeft - slidingCurve.Slices[i].XCenter;
                // Bereken aandrijvend moment : dx * (lamel gewicht - water gewicht)
                moment = moment + horDistance*slidingCurve.Slices[i].Weight;
            }
            return moment;
        }

        /// <summary>
        /// Calculates the Driving Soil Moment on the right side
        /// </summary>
        /// <param name="slidingCurve"></param>
        /// <returns>Driving Soil Moment on the right side</returns>
        private double DrivingSoilMomentRight(SlidingDualCircle slidingCurve)
        {
            double moment = 0;
            for (int i = slidingCurve.RightNumber; i < slidingCurve.Slices.Count; i++)
            {
                double horDistance = xMidRight - slidingCurve.Slices[i].XCenter;
                // Bereken aandrijvend moment : dx * (lamel gewicht - water gewicht)
                moment = moment + horDistance*slidingCurve.Slices[i].Weight;
            }
            return moment;
        }

        /// <summary>
        /// A uplift Van calculation consist of two circel parts and a horicontal part. Equilibrium for the circel parts must include Water against such a part
        /// </summary>
        /// <param name="circle">Circle data</param>
        /// <param name="isLeft">isleft is used to consider if you are calculating on the left or on the right side</param>
        /// <returns>Water moment on circle end</returns>
        private double WaterMomentOnCircleEnd(Circle circle, bool isLeft)
        {
            InterpolationColumn column;
            double zBot = circle.Zcentre - circle.Radius;

            // find the column (within the horizontal part) where the requested data can be found
            if (isLeft)
            {
                column = PreProcessCalculation.GetColumnAtX(circle.Xcentre + Constants.CGeoAccu);
            }
            else
            {
                column = PreProcessCalculation.GetColumnAtX(circle.Xcentre - Constants.CGeoAccu);
            }

            // Find area number in which the bottom of the bar is situated
            int areaNumber = column.DetermineAreaNumber(zBot);
            if (areaNumber < 0)
            {
                return 0.0;
            }

            double poreTop = 0;
            double poreBottom = 0;
            double zTop = 0;
            double zBottom = 0;
            double waterMoment = 0;
            double horizontalInterpolationFactor = (circle.Xcentre - column.XLeft)/(column.XRight - column.XLeft);

            if (areaNumber == 0)
            {
                column.Areas[0].HorizontalInterpolationFactor = horizontalInterpolationFactor;
                zTop = column.Areas[0].ZCoordinateTop();
                zBottom = zBot;
                poreTop = column.Areas[0].PorePressureTop();
                poreBottom = column.Areas[0].PorePressureBottom();
                waterMoment = waterMoment + DetermineWaterMomentOnVertical(poreTop, poreBottom, zTop, zBottom, circle.Zcentre);
            }
            else
            {
                for (int i = 0; i < areaNumber; i++)
                {
                    column.Areas[i].HorizontalInterpolationFactor = horizontalInterpolationFactor;
                    zTop = column.Areas[i].ZCoordinateTop();
                    zBottom = column.Areas[i].ZCoordinateBottom();
                    poreTop = column.Areas[i].PorePressureTop();
                    poreBottom = column.Areas[i].PorePressureBottom();

                    waterMoment = waterMoment + DetermineWaterMomentOnVertical(poreTop, poreBottom, zTop, zBottom, circle.Zcentre);
                }
                // for area number itself
                column.Areas[areaNumber].HorizontalInterpolationFactor = horizontalInterpolationFactor;
                zTop = column.Areas[areaNumber].ZCoordinateTop();
                zBottom = zBot;
                poreTop = column.Areas[areaNumber].PorePressureTop();
                poreBottom = column.Areas[areaNumber].PorePressureAtZ(zBottom);

                waterMoment = waterMoment + DetermineWaterMomentOnVertical(poreTop, poreBottom, zTop, zBottom, circle.Zcentre);
            }
            double result;
            if (isLeft)
            {
                result = -waterMoment;
            }
            else
            {
                result = waterMoment;
            }
            return result;
        }

        /// <summary>
        /// A uplift Van calculation consist of two circel parts and a horicontal part. Equilibrium for the bar parts must include Water against such a part
        /// </summary>
        /// <param name="circle">Circle data: A Bar start or ends at the x coordinate of the circle centre points</param>
        /// <param name="isLeft">isleft is used to consider if you are calculating on the left or on the right sideof a Bar</param>
        /// <returns>The water force at the bar</returns>
        private double WaterForceOnBarEnd(Circle circle, bool isLeft)
        {
            InterpolationColumn column;
            double zBot = circle.Zcentre - circle.Radius;

            // find the column (within the horizontal part) where the requested data can be found
            if (isLeft)
            {
                column = PreProcessCalculation.GetColumnAtX(circle.Xcentre + Constants.CGeoAccu);
            }
            else
            {
                column = PreProcessCalculation.GetColumnAtX(circle.Xcentre - Constants.CGeoAccu);
            }

            // Find area number in which the bottom of the bar is situated
            int areaNumber = column.DetermineAreaNumber(zBot + Constants.CGeoAccu);

            if (areaNumber < 0)
            {
                return 0.0;
            }

            double poreTop = 0;
            double poreBottom = 0;
            double zTop = 0;
            double zBottom = 0;
            double force = 0;
            double horizontalInterpolationFactor = (circle.Xcentre - column.XLeft)/(column.XRight - column.XLeft);

            if (areaNumber == 0)
            {
                column.Areas[0].HorizontalInterpolationFactor = horizontalInterpolationFactor;
                zTop = column.Areas[0].ZCoordinateTop();
                zBottom = zBot;
                poreTop = column.Areas[0].PorePressureTop();
                poreBottom = column.Areas[0].PorePressureBottom();
                force = force + 0.5*(poreTop + poreBottom)*(zTop - zBottom);
            }
            else
            {
                for (int i = 0; i < areaNumber; i++)
                {
                    column.Areas[i].HorizontalInterpolationFactor = horizontalInterpolationFactor;
                    poreTop = column.Areas[i].PorePressureTop();
                    poreBottom = column.Areas[i].PorePressureBottom();
                    zTop = column.Areas[i].ZCoordinateTop();
                    zBottom = column.Areas[i].ZCoordinateBottom();
                    force = force + 0.5*(poreTop + poreBottom)*(zTop - zBottom);
                }

                // for area number itself
                column.Areas[areaNumber].HorizontalInterpolationFactor = horizontalInterpolationFactor;
                zTop = column.Areas[areaNumber].ZCoordinateTop();
                zBottom = zBot;
                poreTop = column.Areas[areaNumber].PorePressureTop();
                poreBottom = column.Areas[areaNumber].PorePressureAtZ(zBottom);
                force = force + 0.5*(poreTop + poreBottom)*(zTop - zBottom);
            }

            return force;
        }

        private double DeterminePoreOnSurface(double aXCoor)
        {
            double result = 0;
            double lTop = PreProcessCalculation.SurfaceLine.Geometry.GetZAtX(aXCoor);
            double yFrea = PreProcessCalculation.Waternet.PhreaticLine.GetZAtX(aXCoor);
            if (yFrea > lTop)
            {
                result = (yFrea - lTop)*PreProcessCalculation.Waternet.UnitWeight;
                result = result*0.5*(yFrea - lTop);
            }
            return result;
        }

        private double DetermineHorizontalWaterForce(SlidingDualCircle slidingCurve)
        {
            var circle = new Circle();
            double lForce = 0.0;
            double lLeftForce = 0.0;
            double lRightForce = 0.0;

            circle.Xcentre = xMidLeft;
            circle.Zcentre = zMidLeft;
            circle.Radius = radLeft;
            lLeftForce = WaterForceOnBarEnd(circle, true);

            circle.Xcentre = xMidRight;
            circle.Zcentre = zMidRight;
            circle.Radius = radRight;
            lRightForce = WaterForceOnBarEnd(circle, false);

            lForce = lLeftForce - lRightForce;
            double leftPoreOnSurface = DeterminePoreOnSurface(xMidLeft);
            Double rightPoreonSurface = DeterminePoreOnSurface(xMidRight);
            return lForce + (leftPoreOnSurface - rightPoreonSurface);
        }

        private double DrivingWaterMomentLeft(SlidingDualCircle slidingCurve)
        {
            double moment = 0;
            for (int i = 0; i <= slidingCurve.LeftNumber; i++)
            {
                // Calculate watermoment due to waterforces on slices
                double horDistance = slidingCurve.Slices[i].XCenter - xMidLeft;
                double vertDistance = zMidLeft - slidingCurve.Slices[i].ZCenterTop;
                moment = moment + horDistance*slidingCurve.Slices[i].VPoreOnSurface +
                         vertDistance*slidingCurve.Slices[i].HPoreOnSurface;
            }
            // Calculate extra moments due to water forces on vertical parts of
            // the geometry within the slip circle area.
            var circle = new Circle();
            {
                circle.Xcentre = slidingCurve.ActiveCircle.X;
                circle.Zcentre = slidingCurve.ActiveCircle.Y;
                circle.Radius = (circle.Zcentre - zTangent);
            }

            moment = moment + WaterMomentOnVerticals(circle);
            // Calculate extra moments due to water forces on vertical parts of
            // the geometry within the slip circle area.
            double waterMoment = WaterMomentOnCircleEnd(circle, true);
            moment = moment + waterMoment;
            return moment;
        }

        private double DrivingWaterMomentRight(SlidingDualCircle slidingCurve)
        {
            double moment = 0;
            for (int i = slidingCurve.RightNumber; i < slidingCurve.Slices.Count; i++)
            {
                // Calculate watermoment due to waterforces on slices
                double horDistance = slidingCurve.Slices[i].XCenter - xMidRight;
                double vertDistance = zMidRight - slidingCurve.Slices[i].ZCenterTop;
                moment = moment + horDistance*slidingCurve.Slices[i].VPoreOnSurface +
                         vertDistance*slidingCurve.Slices[i].HPoreOnSurface;
            }

            // Calculate extra moments due to water forces on vertical parts of
            // the geometry within the slip circle area.

            var circle = new Circle();
            {
                circle.Xcentre = slidingCurve.PassiveCircle.X;
                circle.Zcentre = slidingCurve.PassiveCircle.Y;
                circle.Radius = (circle.Zcentre - zTangent);
            }

            moment = moment + WaterMomentOnVerticals(circle);
            moment = moment + WaterMomentOnCircleEnd(circle, false);

            return moment;
        }

        private double UniformLoadMoment(double Xcentre, double Xstart, double Xend)
        {
            double loadMoment = 0;
            foreach (var uniformLoad in PreProcessCalculation.UniformLoads)
            {
                loadMoment += uniformLoad.GetMomentRoundX(Xcentre, Xstart, Xend);
            }
            return loadMoment;
        }

        private void DetermineDrivingMoments(SlidingDualCircle slidingCurve)
        {
            totDrivingLeft = 0.0;
            totDrivingRight = 0.0;
            // calculate driving moments
            soilMomentLeft = DrivingSoilMomentLeft(slidingCurve);
            soilMomentRight = DrivingSoilMomentRight(slidingCurve);
            // watermoments from free water
            watermomentLeft = DrivingWaterMomentLeft(slidingCurve);
            watermomentRight = DrivingWaterMomentRight(slidingCurve);
            // load moments if present
            loadMomentLeft = UniformLoadMoment(xMidLeft, slidingCurve.LeftPoint.X, xMidLeft);
            loadMomentRight = UniformLoadMoment(xMidRight, xMidRight, slidingCurve.RightPoint.X);
            // the total driving moments
            totDrivingLeft = soilMomentLeft + watermomentLeft + loadMomentLeft;
            totDrivingRight = soilMomentRight + watermomentRight + loadMomentRight;
            reverseGeom = Math.Abs(totDrivingRight) > Math.Abs(totDrivingLeft);
        }

        /// <summary>
        /// Calculates the horizontal shear force between two circle parts
        /// </summary>
        /// <param name="slidingCurve"></param>
        /// <param name="startSlice"></param>
        /// <param name="endSlice"></param>
        /// <param name="factorFo"></param>
        /// <returns></returns>
        private double GetHorResistingForce(SlidingDualCircle slidingCurve, int startSlice, int endSlice, double factorFo)
        {
            const double pir = Math.PI/180.0;
            double som = 0;
//            double sinPhi = 0;
//            double cosPhi = 0;
            double tau = 0;

            for (int i = startSlice; i <= endSlice; i++)
            {
                Slice slice = slidingCurve.Slices[i];
                double coh = slice.Cohesion;

                // cohesie
                //sinPhi = Math.Sin(slice.Phi*pir);
                //cosPhi = Math.Cos(slice.Phi*pir);
                //LCosDilet = Math.Cos(LSlice.Dilatancy * Pir);
                //LSinDilet = Math.Sin(LSlice.Dilatancy * Pir);
                //LSigTerm = LCosDilet * LSinPhi / (1 - LSinDilet * LSinPhi);
                //LCohTerm = LCosDilet * LCosPhi / (1 - LSinDilet * LSinPhi);
                // bereken de coh en phi uit de spannings  tabel

                //LTau = coh * LCohTerm + LSlice.EffectiveStress * LSigTerm;
                tau = coh + slice.EffectiveStress*Math.Tan(slice.Phi*pir);
                som = som + slice.ArcLength*tau;
                // Nieuwe Tau - alfa
                slice.NormalStress = slice.EffectiveStress;
                slice.ShearStress = tau/factorFo;
            }
            return som;
        }

        /// <summary>
        /// Calculates the resisting moments of a circle part
        /// </summary>
        /// <param name="slidingCurve"></param>
        /// <param name="startSlice"></param>
        /// <param name="endSlice"></param>
        /// <param name="radius"></param>
        /// <param name="factorFo"></param>
        /// <returns></returns>
        private double GetResistingMoment(SlidingDualCircle slidingCurve, int startSlice, int endSlice, double radius, double factorFo)
        {
            const double Pir = Math.PI/180.0;
            int i;
            double coh;
            double tgPhi;
            double Salfa;
            double alfa;
            double h;
            double tgalfa;
            double SomMoment;
            //           double LCohTerm;
/*
            double LSigTerm;
            double LSinPhi;
            double LCosPhi;
            double LCosDilet;
            double LSinDilet;
*/
            double LTau;
            SomMoment = 0.0;
            for (i = startSlice; i <= endSlice; i++)
            {
                coh = slidingCurve.Slices[i].Cohesion;
                // cohesie
/*
                LSinPhi = Math.Sin(slidingCurve.Slices[i].Phi * Pir);
                LCosPhi = Math.Cos(slidingCurve.Slices[i].Phi * Pir);
                LCosDilet = Math.Cos(slidingCurve.Slices[i].Dilatancy * Pir);
                LSinDilet = Math.Sin(slidingCurve.Slices[i].Dilatancy * Pir);
*/
                tgPhi = Math.Tan(slidingCurve.Slices[i].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);
*/
                alfa = slidingCurve.Slices[i].BottomAngle;
                // Hoek onderkant lamel [rad]
                if (DrivingMoment > 0.0)
                {
                    alfa = -alfa;
                }
                // Passieve alfa nagatief maken
                // Afsnuiten indien nodig
                h = 0.5*Math.Atan(tgPhi/factorFo) - 0.25*Math.PI;
                if (alfa <= h)
                {
                    alfa = h;
                }
                tgalfa = Math.Tan(alfa);
                // Bereken het tegen werkend  moment
                LTau = coh + slidingCurve.Slices[i].EffectiveStress*tgPhi;
                SomMoment = SomMoment + slidingCurve.Slices[i].ArcLength*LTau/(1 + tgalfa*tgPhi/factorFo);
                Salfa = (factorFo*slidingCurve.Slices[i].EffectiveStress - coh*tgalfa)/(factorFo + tgPhi*tgalfa);
                slidingCurve.Slices[i].NormalStress = Salfa;
                slidingCurve.Slices[i].ShearStress = (coh + Salfa*tgPhi)/factorFo;
            }
            return SomMoment*radius;
        }

        /// <summary>
        /// Calculates the safety factor for Uplift-Van method
        /// </summary>
        /// <param name="slidingCurve">The considered sliding curve (2 circles + horizontal part)</param>
        /// <returns>Safety factor for Uplift-Van method</returns>
        private double LiftMethode(SlidingDualCircle slidingCurve)
        {
            // Initialisatie
            resistingMom = 0.0;
            double fn = StartValueSafetyfactor;
            bool ready = false;
            var numberOfIterations = 0;

            double armLeft = radLeft - 1d/3d*(slidingCurve.Slices[slidingCurve.LeftNumber].TopRight.Y - slidingCurve.Slices[slidingCurve.LeftNumber].BottomRight.Y);
            double armRight = radRight - 1d/3d*(slidingCurve.Slices[slidingCurve.RightNumber].TopLeft.Y - slidingCurve.Slices[slidingCurve.RightNumber].BottomLeft.Y);

            drivingMoment = totDrivingLeft/armLeft + totDrivingRight/armRight;

            double horWaterForce = DetermineHorizontalWaterForce(slidingCurve);
            // Driving moment is actualy a force 
            drivingMoment = drivingMoment + horWaterForce;

            // Check driving moment on 1 kNm
            if ((Math.Abs(drivingMoment) < 0.001))
            {
                ready = true;
            }
            // Start iteratie
            if (!ready)
            {
                do
                {
                    numberOfIterations++;
                    double fo = fn;
                    resistingMom = 0.0;

                    // bereken de momenten van de linker en rechter helfr
                    // rekening houdend met de draairichting (reversegeom of niet )
                    double leftResMom = GetResistingMoment(slidingCurve, 0, slidingCurve.LeftNumber, radLeft, fo);
                    double rightResMom = GetResistingMoment(slidingCurve, slidingCurve.RightNumber, slidingCurve.Slices.Count - 1, radRight, fo);
                    double horResForce = GetHorResistingForce(slidingCurve, slidingCurve.LeftNumber + 1, slidingCurve.RightNumber - 1, fo);

                    resistingMom = Math.Abs(leftResMom)/armLeft + Math.Abs(rightResMom)/armRight + Math.Abs(horResForce);
                    // bereken de actieve passieve en hor kracht
                    if (reverseGeom)
                    {
                        passiveForce = (Math.Abs(totDrivingLeft) + Math.Abs(leftResMom/fo))/armLeft;
                        activeForce = (Math.Abs(totDrivingRight) - Math.Abs(rightResMom/fo))/armRight;
                        activeDrivingDMoment = totDrivingRight;
                        activeResistingMoment = rightResMom;
                        passiveDrivingMoment = totDrivingLeft;
                        passiveResistingMoment = leftResMom;
                    }
                    else
                    {
                        activeForce = (Math.Abs(totDrivingLeft) - Math.Abs(leftResMom)/fo)/armLeft;
                        passiveForce = (Math.Abs(totDrivingRight) + Math.Abs(rightResMom/fo))/armRight;
                        activeDrivingMoment = totDrivingLeft;
                        activeResistingMoment = leftResMom;
                        passiveDrivingMoment = totDrivingRight;
                        passiveResistingMoment = rightResMom;
                    }

                    horForce = horResForce/fo;
                    if ((numberOfIterations == 1))
                    {
                        slidingCurve.ActiveForce0 = activeForce;
                        slidingCurve.PassiveForce0 = passiveForce;
                        slidingCurve.HorizontalForce0 = horForce;
                        resistingMom0 = resistingMom;
                    }
                    else
                    {
                        slidingCurve.ActiveForce = activeForce;
                        slidingCurve.PassiveForce = passiveForce;
                        slidingCurve.HorizontalForce = horForce;
                        slidingCurve.DrivingMomentActive = activeDrivingMoment;
                        slidingCurve.DrivingMomentPassive = passiveDrivingMoment;
                        slidingCurve.ResistingMomentActive = activeResistingMoment;
                        slidingCurve.ResistingMomentPassive = passiveResistingMoment;
                    }

                    fn = resistingMom/Math.Abs(drivingMoment);

                    // 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(resistingMom) < 0.001) || (fn < 0.001))
                    {
                        ready = true;
                    }
                } while (!(ready));
            }

            double result = fn;
            if (Math.Abs(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 " + Convert.ToString(Constants.CMaxIter) + " iterations"));
                }
                else if ((Math.Abs(resistingMom) < 0.001) || (fn < 0.001))
                {
                    result = 1001;
                    LogManager.Add(new LogMessage(LogMessageType.Warning, this, "Resisting moment too small"));
                }
            }
            return result;
        }

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

        private bool IfMinimalPlanedepthIsNotOk(SlidingDualCircle slidingCurve)
        {
            for (int i = 0; i < slidingCurve.Slices.Count; i++)
            {
                if ((slidingCurve.Slices[i].ZCenterTop - slidingCurve.Slices[i].ZCenterBottom) >
                    PreProcessCalculation.Constraints.SlipPlaneMinDepth)
                {
                    return false;
                }
            }
            return true;
        }

        private bool IfMinimalPlaneLengthIsNotOk(SlidingDualCircle slidingCurve)
        {
            double som = 0;
            for (int i = 0; i < slidingCurve.Slices.Count; i++)
            {
                som = som + slidingCurve.Slices[i].ArcLength;
            }
            if (som > PreProcessCalculation.Constraints.SlipPlaneMinLength)
            {
                return false;
            }
            return true;
        }

        private double GetStabilityFactor(SlidingDualCircle slidingCurve)
        {
            // Calculate the driving moments
            DetermineDrivingMoments(slidingCurve);

            // Calculate stability factor
            slidingCurve.SafetyFactor = LiftMethode(slidingCurve);

            return slidingCurve.SafetyFactor;
        }

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

        private double GetActiveRadius()
        {
            double result;
            result = Math.Abs(ZActiveCentre - ZTangent);
            return result;
        }

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

        private double GetPassiveRadius()
        {
            double result;
            result = Math.Abs(ZPassiveCentre - ZTangent);
            return result;
        }

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

        //===========================================================================
    }
}