using System;
using System.Collections.Generic;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Mathematics;
using Deltares.Standard;
using Deltares.Standard.Extensions;

namespace Deltares.Stability.Calculation2
{
    public class SlicesCreator
    {
        private List<GeometryPoint> changingPoints;

        private SoilProfile2D profile2D;

        // point of plane to be calculated
        private GeometryPointString spencerplane;
        private double striveWidth; // gewenste breedte van de slices
        private GeometryPointString surfaceLine;
        private Waternet waternet;

        public GeometryPointString Spencerplane
        {
            get
            {
                return spencerplane;
            }
            set
            {
                spencerplane = value;
            }
        }

        public double StriveWidth
        {
            get
            {
                return striveWidth;
            }
            set
            {
                striveWidth = value;
            }
        }

        public List<GeometryPoint> ChangingPoints
        {
            get
            {
                return changingPoints;
            }
            set
            {
                changingPoints = value;
            }
        }

        public void SetSoilProfile2D(SoilProfile2D profile)
        {
            profile2D = profile;
        }

        public void SetWaternet(Waternet waternetExt)
        {
            waternet = waternetExt;
        }

        public void SetSurfaceLine(SurfaceLine2 line)
        {
            surfaceLine = line.Geometry;
        }

        public static double Distance(double x1, double z1, double x2, double z2)
        {
            return Math.Sqrt((x1 - x2)*(x1 - x2) + (z1 - z2)*(z1 - z2));
        }

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

        public double DetermineSpencerBottomPoint(double AXCoor)
        {
            double result;
            ;
            bool Found;
            int i;
            i = 0;
            result = -999.99;
            do
            {
                i++;
                Found = ((AXCoor >= spencerplane.Points[i - 1].X) && (AXCoor <= spencerplane.Points[i].X));
                if (Found)
                {
                    result = MStabDatafunctions.LinInpolY(spencerplane.Points[i - 1].X, spencerplane.Points[i - 1].Z,
                                                          spencerplane.Points[i].X, spencerplane.Points[i].Z, AXCoor);
                }
            } while (!(Found || (i == spencerplane.Points.Count - 1)));

            return result;
        }

        public static bool aboveLine(GeometryPointString APlane, GeometryPoint APoint)
        {
            bool result = false;
            List<double> zCoors = APlane.GetAllZatXForSurface(APoint.X);
            if (zCoors.Count > 0)
            {
                for (int i = 0; i < zCoors.Count; i++)
                {
                    if (zCoors[i] <= APoint.Z)
                    {
                        result = true;
                        return result;
                    }
                }
            }
            return result;
        }

        public void AddPlaneIntersectionPointsWithLayer(SlidingCurve slidingCurve)
        {
            for (int i = 0; i < profile2D.Surfaces.Count; i++)
            {
                List<GeometryPoint> intersections =
                    spencerplane.IntersectionXzPointsWithGeometryString(profile2D.Surfaces[i].GeometrySurface.OuterLoop);
                {
                    for (int j = 0; j < intersections.Count; j++)
                    {
                        NextSlice(slidingCurve, intersections[j].X);
                    }
                }
            }
        }

        public void AddPlaneIntersectionPointsWithFreaticLine(SlidingCurve slidingCurve)
        {
            if (waternet.PhreaticLine != null)
            {
                List<GeometryPoint> intersections =
                    spencerplane.IntersectionXzPointsWithGeometryString(waternet.PhreaticLine);
                for (int j = 0; j < intersections.Count; j++)
                {
                    NextSlice(slidingCurve, intersections[j].X);
                }
            }
        }

        public void AddPlaneIntersectionPointsWithWaternetLinesFreaticLine(SlidingCurve slidingCurve)
        {
            foreach (var waternetLine in waternet.WaternetLineList)

            {
                List<GeometryPoint> intersections = spencerplane.IntersectionXzPointsWithGeometryString(waternetLine);
                for (int j = 0; j < intersections.Count; j++)
                {
                    NextSlice(slidingCurve, intersections[j].X);
                }
            }
        }

        public void CreateSpencerSlices(SlidingCurve slidingCurve)
        {
            // first the points from the slipplane
            slidingCurve.LeftPoint = spencerplane.Points[0];
            slidingCurve.RightPoint = spencerplane.Points[spencerplane.Points.Count - 1];
            for (int i = 0; i < spencerplane.Points.Count; i++)
            {
                NextSlice(slidingCurve, spencerplane.Points[i].X);
            }

            for (int i = 0; i < changingPoints.Count; i++)
            {
                // Only relevant geomery points
                if ((changingPoints[i].X > slidingCurve.LeftPoint.X) && (changingPoints[i].X < slidingCurve.RightPoint.X))
                {
                    if (aboveLine(spencerplane, changingPoints[i]))
                    {
                        NextSlice(slidingCurve, changingPoints[i].X);
                    }
                }
            }

            // Intersection points with layer boundaries
            AddPlaneIntersectionPointsWithLayer(slidingCurve);
            // Intersection with freatic line
            if (waternet.PhreaticLine != null)
            {
                AddPlaneIntersectionPointsWithFreaticLine(slidingCurve);
            }
            // intersection with water net lines
            if (waternet.WaternetLineList.Count > 0)
            {
                AddPlaneIntersectionPointsWithWaternetLinesFreaticLine(slidingCurve);
            }

            // sorteer de x koordinaten van de lamellen
            slidingCurve.Slices.Sort();
            // Divide large slices in smaller ones
            DivideLargeSlicesInSmaller(slidingCurve);
            // Sort slices again
            slidingCurve.Slices.Sort();
            // Fill the Right side too
            FillRightSide(slidingCurve);
            // gooi de laatste dumy lamel weg
            slidingCurve.Slices.RemoveAt(slidingCurve.Slices.Count - 1);
        }

        public void CreateBishopSlices(SlidingCircle slidingCircle)
        {
            if (DetermineSurfacePoints(slidingCircle))
            {
                NextSlice(slidingCircle, slidingCircle.LeftPoint.X);
                NextSlice(slidingCircle, slidingCircle.RightPoint.X);
                // Add the centre of the circle
                NextSlice(slidingCircle, slidingCircle.CircleCenterX);
                bool outOfGeometry = false;
                if (surfaceLine.Points.Count > 1)
                {
                    outOfGeometry = (slidingCircle.LeftPoint.X < surfaceLine.Points[0].X) ||
                                    (slidingCircle.RightPoint.X >
                                     surfaceLine.Points[surfaceLine.Points.Count - 1].X);
                }
                // Punten uit de puntentabel (die laagscheidingen/PN-lijnen vormen
                // Deze punten zijn bij het preprocessen gemaakt
                double radius = slidingCircle.Radius;

                for (int i = 0; i < changingPoints.Count; i++)
                {
                    // Only relevant geomery points
                    if ((changingPoints[i].X > slidingCircle.LeftPoint.X) && (changingPoints[i].X < slidingCircle.RightPoint.X))
                    {
                        if ((Distance(changingPoints[i].X, changingPoints[i].Z, slidingCircle.Center.X, slidingCircle.Center.Y) <= radius + Constants.CGeoAccu))
                        {
                            NextSlice(slidingCircle, changingPoints[i].X);
                        }
                    }
                }
                // Intersection points with layer boundaries
                DetermineCircleIntersectionPointsWithLayers(slidingCircle);
                // Intersection with freatic line
                if (waternet.PhreaticLine != null)
                {
                    DetermineCircleIntersectionPointsWithFreaticLine(slidingCircle);
                }
                // Intersection with waternet lines
                if (waternet.WaternetLineList.Count > 0)
                {
                    DetermineCircleIntersectionPointsWithWaternetLines(slidingCircle);
                }

                if (outOfGeometry)
                {
                    DetermineCircleIntersectionPointOutOfGeometry(slidingCircle);
                }

                // sorteer de x koordinaten van de lamellen
                slidingCircle.Slices.Sort();
                // Divide large slices in smaller ones
                DivideLargeSlicesInSmaller(slidingCircle);
                // Sort slices again
                slidingCircle.Slices.Sort();
                // Fill the Right side too
                FillRightSide(slidingCircle);
                // gooi de laatste dumy lamel weg
                if (slidingCircle.Slices.Count > 0)
                {
                    slidingCircle.Slices.RemoveAt(slidingCircle.Slices.Count - 1);
                }
                // calculate additional geometry data  CALCULATED ON THE FLY
                //for (int i = 0; i < slices.Count; i++)
                //{
                //    slices[i].CreateGeometryItems();
                //}
            }
        }

        /// <summary>
        /// First find the points at left and right of the slipplane
        /// Next se what pre defined points are in the area of the slipplane
        /// finally find the intersection points with layers and freatic line
        /// </summary>
        /// <param name="slidingCircle"></param>
        public void CreateUpliftVanSlices(SlidingDualCircle slidingCircle)
        {
            if (DetermineSurfacePoints(slidingCircle))
            {
                addSurfacePoints(slidingCircle);

                addCentrePoints(slidingCircle);

                // Punten uit de puntentabel (die laagscheidingen/PN-lijnen vormen
                addPoinstAboveVanSlipPlane(slidingCircle);

                DetermineUpliftIntersectionPointsWithLayers(slidingCircle);
                DetermineUpliftIntersectionPointsWithWaternetlines(slidingCircle);

                // met headlines is eigenlijk niet nodig maar gebeurt om te vergelijken met oude versie
                // DetermineUpliftIntersectionPointsWithHeadlines(slidingCircle);

                // sorteer de x koordinaten van de lamellen
                slidingCircle.Slices.Sort();
                // Divide large slices in smaller ones
                DivideLargeSlicesInSmaller(slidingCircle);
                // Sort slices again
                slidingCircle.Slices.Sort();
                // Fill the Right side too
                FillRightSide(slidingCircle);
                // gooi de laatste dumy lamel weg
                if (slidingCircle.Slices.Count > 0)
                {
                    slidingCircle.Slices.RemoveAt(slidingCircle.Slices.Count - 1);
                }
                // caluclat aditional geometry data  CALCULATED ON THE FLY
            }
        }

        private double DetermineUpliftBottompoint(SlidingDualCircle slidingCurve, double x)
        {
            Point2D leftCircle = slidingCurve.ActiveCircle.X > slidingCurve.PassiveCircle.X
                                     ? slidingCurve.PassiveCircle
                                     : slidingCurve.ActiveCircle;
            Point2D rightCircle = slidingCurve.ActiveCircle.X > slidingCurve.PassiveCircle.X
                                      ? slidingCurve.ActiveCircle
                                      : slidingCurve.PassiveCircle;

            if (slidingCurve.HasWig == false)
                //if (stabilityModel.ModelOption != ModelOptions.UpliftSpencer)
            {
                if (x < leftCircle.X)
                {
                    return MStabDatafunctions.ZAtXOnCircle(x, leftCircle.X, leftCircle.Y,
                                                           leftCircle.Y - slidingCurve.TangentLine);
                }
                else
                {
                    if (x > rightCircle.X)
                    {
                        return MStabDatafunctions.ZAtXOnCircle(x, rightCircle.X, rightCircle.Y,
                                                               rightCircle.Y - slidingCurve.TangentLine);
                    }
                    else
                    {
                        return slidingCurve.TangentLine;
                    }
                }
            }
            else
            {
                // Lift spencer has a wig at passive side
                if (x < leftCircle.X)
                {
                    if (leftCircle == slidingCurve.ActiveCircle)
                    {
                        return MStabDatafunctions.ZAtXOnCircle(x, leftCircle.X, leftCircle.Y,
                                                               leftCircle.Y - slidingCurve.TangentLine);
                    }
                    else
                    {
                        double zSurface = MStabDatafunctions.GetZTopOrBottomAtSurface(slidingCurve.LeftPoint.X,
                                                                                      surfaceLine, true);
                        return MStabDatafunctions.LinInpolY(slidingCurve.LeftPoint.X, zSurface, leftCircle.X,
                                                            slidingCurve.TangentLine, x);
                    }
                }
                else
                {
                    if (x > rightCircle.X)
                    {
                        if (leftCircle == slidingCurve.ActiveCircle)
                        {
                            double zSurface = MStabDatafunctions.GetZTopOrBottomAtSurface(slidingCurve.RightPoint.X,
                                                                                          surfaceLine, true);
                            return MStabDatafunctions.LinInpolY(slidingCurve.RightPoint.X, zSurface, rightCircle.X,
                                                                slidingCurve.TangentLine, x);
                        }
                        else
                        {
                            return MStabDatafunctions.ZAtXOnCircle(x, rightCircle.X, rightCircle.Y,
                                                                   rightCircle.Y - slidingCurve.TangentLine);
                        }
                    }
                    else
                    {
                        return slidingCurve.TangentLine;
                    }
                }
            }
        }

        private double DetermineZBottom(SlidingCurve slidingCurve, double x)
        {
            if (slidingCurve is SlidingCircle)
            {
                var slidingCircle = (SlidingCircle) slidingCurve;
                return MStabDatafunctions.ZAtXOnCircle(x, slidingCircle.Center.X, slidingCircle.Center.Y,
                                                       slidingCircle.Radius);
            }
            else if (slidingCurve is SlidingDualCircle)
            {
                var slidingDualCircle = (SlidingDualCircle) slidingCurve;
                return DetermineUpliftBottompoint(slidingDualCircle, x);
            }
            else if (slidingCurve is SlidingPlane)
            {
                return DetermineSpencerBottomPoint(x);
            }
            else
            {
                throw new NotImplementedException("DetermineZBottom for " + slidingCurve.GetType());
            }
        }

        private void FillRightSide(SlidingCurve slidingCurve)
        {
            int i;
            for (i = 0; i < slidingCurve.Slices.Count - 1; i++)
            {
                slidingCurve.Slices[i].TopRight.X = slidingCurve.Slices[i + 1].TopLeft.X;
                slidingCurve.Slices[i].TopRight.Y = slidingCurve.Slices[i + 1].TopLeft.Y;
                slidingCurve.Slices[i].BottomRight.X = slidingCurve.Slices[i + 1].BottomLeft.X;
                slidingCurve.Slices[i].BottomRight.Y = slidingCurve.Slices[i + 1].BottomLeft.Y;
            }
        }

        /// <summary>
        /// Add a slice if betwen the boundaries 
        /// and not les than CMinSliceWidth from an existing slice
        /// </summary>
        /// <param name="slidingCurve"></param>
        /// <param name="x"></param>
        private void NextSlice(SlidingCurve slidingCurve, double x)
        {
            // First check if x is in slipplane area.
            if ((x > slidingCurve.LeftPoint.X - Slice.CMinSliceWidth) && (x < slidingCurve.RightPoint.X + Slice.CMinSliceWidth))
            {
                // Kijk of het een nieuwe X-koordinaat betreft
                bool found = false;
                foreach (var existingSlice in slidingCurve.Slices)
                {
                    if (Math.Abs(x - existingSlice.TopLeft.X) < Slice.CMinSliceWidth)
                    {
                        found = true;
                    }
                }

                // Add X if not yet present
                if (!found)
                {
                    var slice = new Slice();
                    slice.TopLeft.X = x;
                    slice.TopLeft.Y = surfaceLine.GetZAtX(x);
                    slice.BottomLeft.X = x;
                    slice.BottomLeft.Y = DetermineZBottom(slidingCurve, x);
                    slidingCurve.Slices.Add(slice);
                }
            }
        }

        /// <summary>
        /// This procedure finds the valid slip circle is case that more than 2 intersections points between the slip circle and the surface line were found.
        ///  </summary>
        /// <param name="slidingCircle">The current slip circle</param>
        /// <param name="intersectionPoints">The list of intersections points between the slip circle and the surface line</param>
        /// <returns>True means that a valid slip circle was found</returns>
        public bool HandleIntersectionsMoreThanTwo(SlidingCircle slidingCircle, IList<GeometryPoint> intersectionPoints)
        {
            double dx = 0;
            bool result = false;
            // The intersection points situated above the centre of the slip circle are not used because leading to unrealistic slip plane
            var validIntersectionPoints = new List<GeometryPoint>();
            for (int i = 0; i < intersectionPoints.Count; i++)
            {
                if (intersectionPoints[i].Z <= slidingCircle.Center.Y)
                {
                    validIntersectionPoints.Add(intersectionPoints[i]);
                }
            }
            for (int i = 1; i < validIntersectionPoints.Count; i++)
            {
                double xmid = 0.5 * (validIntersectionPoints[i].X + validIntersectionPoints[i - 1].X);
                double zmid = surfaceLine.GetZAtX(xmid);
                // A valid slip circle is  situated completely below the surface line
                if (Distance(xmid, zmid, slidingCircle.Center.X, slidingCircle.Center.Y) < slidingCircle.Radius)
                {
                    if (Math.Abs(validIntersectionPoints[i].X - validIntersectionPoints[i - 1].X) > dx)
                    {
                        slidingCircle.LeftPoint = new GeometryPoint(validIntersectionPoints[i - 1].X, 0, validIntersectionPoints[i - 1].Z);
                        slidingCircle.RightPoint = new GeometryPoint(validIntersectionPoints[i].X, 0, validIntersectionPoints[i].Z);
                        dx = Math.Abs(validIntersectionPoints[i].X - validIntersectionPoints[i - 1].X);
                        result = true;
                    }
                }
            }

            return result;
        }

        private void DetermineCircleIntersectionPointOutOfGeometry(SlidingCircle slidingCircle)
        {
            bool leftside = (slidingCircle.LeftPoint.X < surfaceLine.Points[0].X);
            SoilProfile1D soil1D;
            double bounX;
            if (leftside)
            {
                bounX = surfaceLine.Points[0].X;
                soil1D = profile2D.GetSoilProfile1D(bounX);
                for (int i = 0; i < soil1D.Layers.Count; i++)
                {
                    //   NextSlice(slidingCircle, interpoints[i].X);
                }
            }
            else
            {
                bounX = surfaceLine.Points[surfaceLine.Points.Count - 1].X;
                soil1D = profile2D.GetSoilProfile1D(bounX);
                for (int i = 0; i < soil1D.Layers.Count; i++) {}
            }
        }

        private void DetermineCircleIntersectionPointsWithLayers(SlidingCircle slidingCircle)
        {
            foreach (var surface in profile2D.Surfaces)
            {
                bool intersectionIsNotPossible;
                intersectionIsNotPossible = (surface.GeometrySurface.OuterLoop.GetMaxX() < slidingCircle.LeftPoint.X) ||
                                            (surface.GeometrySurface.OuterLoop.GetMinX() > slidingCircle.RightPoint.X) ||
                                            (surface.GeometrySurface.OuterLoop.GetMaxZ() <
                                             slidingCircle.Center.Y - slidingCircle.Radius);

                if (!intersectionIsNotPossible)
                {
                    IList<GeometryPoint> interpoints =
                        LineHelper.IntersectionPointsWithCircle(slidingCircle.Center.X,
                                                                slidingCircle.Center.Y, slidingCircle.Radius, surface.GeometrySurface.OuterLoop.Points);
                    for (int i = 0; i < interpoints.Count; i++)
                    {
                        if ((interpoints[i].X > slidingCircle.LeftPoint.X) &&
                            (interpoints[i].X < slidingCircle.RightPoint.X) &&
                            (interpoints[i].Z > slidingCircle.Center.Y - slidingCircle.Radius))
                        {
                            NextSlice(slidingCircle, interpoints[i].X);
                        }
                    }
                }
            }
        }

        private void DetermineCircleIntersectionPointsWithFreaticLine(SlidingCircle slidingCircle)
        {
            if (waternet.PhreaticLine != null)
            {
                IList<GeometryPoint> interpoints =
                    LineHelper.IntersectionPointsWithCircle(slidingCircle.Center.X,
                                                            slidingCircle.Center.Y, slidingCircle.Radius, waternet.PhreaticLine.Points);
                for (int i = 0; i < interpoints.Count; i++)
                {
                    if ((interpoints[i].X > slidingCircle.LeftPoint.X) &&
                        (interpoints[i].X < slidingCircle.RightPoint.X))
                    {
                        NextSlice(slidingCircle, interpoints[i].X);
                    }
                }
            }
        }

        private void DetermineCircleIntersectionPointsWithWaternetLines(SlidingCircle slidingCircle)
        {
            foreach (var waternetLine in waternet.WaternetLineList)

            {
                IList<GeometryPoint> interpoints = LineHelper.IntersectionPointsWithCircle(
                    slidingCircle.Center.X, slidingCircle.Center.Y, slidingCircle.Radius, waternetLine.Points);
                for (int i = 0; i < interpoints.Count; i++)
                {
                    if ((interpoints[i].X > slidingCircle.LeftPoint.X) &&
                        (interpoints[i].X < slidingCircle.RightPoint.X))
                    {
                        NextSlice(slidingCircle, interpoints[i].X);
                    }
                }
            }
        }

        /// <summary>
        /// Determine the surface points for the uplift van calculation. 
        /// </summary>
        /// <param name="upliftVanSlidingCurve">The current Uplift Van sliding curve</param>
        /// <returns>returns true if the left and right surface point is added to the sliding curve</returns>
        private bool DetermineSurfacePoints(SlidingDualCircle upliftVanSlidingCurve)
        {
            bool result;
            var activeIntersectionPoints = LineHelper.ExtendedSurfaceIntersectionPointsWithCircle
                (upliftVanSlidingCurve.ActiveCircle.X, upliftVanSlidingCurve.ActiveCircle.Y,
                 upliftVanSlidingCurve.ActiveCircle.Y - upliftVanSlidingCurve.TangentLine, surfaceLine.Points);
            var passiveIntersectionPoints = LineHelper.ExtendedSurfaceIntersectionPointsWithCircle
                (upliftVanSlidingCurve.PassiveCircle.X, upliftVanSlidingCurve.PassiveCircle.Y,
                 upliftVanSlidingCurve.PassiveCircle.Y - upliftVanSlidingCurve.TangentLine, surfaceLine.Points);
            var horizontalPart = new Line();
            horizontalPart.BeginPoint = new GeometryPoint(upliftVanSlidingCurve.ActiveCircleCenterX, 0, upliftVanSlidingCurve.TangentLine);
            horizontalPart.EndPoint = new GeometryPoint(upliftVanSlidingCurve.PassiveCircleCenterX, 0, upliftVanSlidingCurve.TangentLine);
            var horizontalIntersectionPoints = surfaceLine.IntersectionPointsXzWithLineXzWithAllPoints(horizontalPart);
            horizontalIntersectionPoints.Sort();

            if ((activeIntersectionPoints.Count == 0) || (passiveIntersectionPoints.Count == 0))
            {
                return false;
            }

            result = false;
            double minimumXDistance = double.MaxValue;
            var currentLeftPoint = new GeometryPoint();
            var currentRightPoint = new GeometryPoint();
            // Case 1: the Active circle is at left side
            if (upliftVanSlidingCurve.ActiveCircle.X <= upliftVanSlidingCurve.PassiveCircle.X + Constants.CSameSliceDist)
            {
                for (int i = 0; i < activeIntersectionPoints.Count; i++)
                {
                    // Search the intersection point at the left side of the active circle
                    if (activeIntersectionPoints[i].X < upliftVanSlidingCurve.ActiveCircle.X)
                    {
                        // If several intersection points are found, use the closest to the Xcenter 
                        if (Math.Abs(activeIntersectionPoints[i].X - upliftVanSlidingCurve.ActiveCircle.X) <
                            minimumXDistance)
                        {
                            currentLeftPoint = activeIntersectionPoints[i];
                            minimumXDistance =
                                Math.Abs(activeIntersectionPoints[i].X - upliftVanSlidingCurve.ActiveCircle.X);
                        }
                        result = true;
                    }
                }
                upliftVanSlidingCurve.LeftPoint = new GeometryPoint(currentLeftPoint.X, 0, currentLeftPoint.Z);

                minimumXDistance = double.MaxValue;
                if (result)
                {
                    result = false;
                    // The RightPoint can be situated in the horizontal part, this is not allowed
                    if (horizontalIntersectionPoints.Count == 0)
                    {
                        for (int i = 0; i < passiveIntersectionPoints.Count; i++)
                        {
                            if (passiveIntersectionPoints[i].X > upliftVanSlidingCurve.PassiveCircle.X)
                            {
                                // If several intersection points are found, use the closest to the Xcenter 
                                if (Math.Abs(passiveIntersectionPoints[i].X - upliftVanSlidingCurve.PassiveCircle.X) < minimumXDistance)
                                {
                                    currentRightPoint = passiveIntersectionPoints[i];
                                    minimumXDistance = Math.Abs(passiveIntersectionPoints[i].X - upliftVanSlidingCurve.PassiveCircle.X);
                                }
                                result = true;
                            }
                        }
                        upliftVanSlidingCurve.RightPoint = new GeometryPoint(currentRightPoint.X, 0, currentRightPoint.Z);
                    }
                }
            }
                // Case 2: the Active circle is at right side
            else
            {
                for (int i = 0; i < passiveIntersectionPoints.Count; i++)
                {
                    if (passiveIntersectionPoints[i].X < upliftVanSlidingCurve.PassiveCircle.X)
                    {
                        // If several intersection points are found, use the closest to the Xcenter 
                        if (Math.Abs(passiveIntersectionPoints[i].X - upliftVanSlidingCurve.PassiveCircle.X) < minimumXDistance)
                        {
                            currentLeftPoint = passiveIntersectionPoints[i];
                            minimumXDistance = Math.Abs(passiveIntersectionPoints[i].X - upliftVanSlidingCurve.PassiveCircle.X);
                        }
                        result = true;
                    }
                }
                upliftVanSlidingCurve.LeftPoint = new GeometryPoint(currentLeftPoint.X, 0, currentLeftPoint.Z);

                minimumXDistance = double.MaxValue;
                if (result)
                {
                    result = false;
                    // The LeftPoint can be situated in the horizontal part, this is not allowed
                    if (horizontalIntersectionPoints.Count == 0)
                    {
                        for (int i = 0; i < activeIntersectionPoints.Count; i++)
                        {
                            if (activeIntersectionPoints[i].X > upliftVanSlidingCurve.ActiveCircle.X)
                            {
                                // If several intersection points are found, use the closest to the Xcenter 
                                if (Math.Abs(activeIntersectionPoints[i].X - upliftVanSlidingCurve.PassiveCircle.X) < minimumXDistance)
                                {
                                    currentRightPoint = activeIntersectionPoints[i];
                                    minimumXDistance = Math.Abs(activeIntersectionPoints[i].X - upliftVanSlidingCurve.PassiveCircle.X);
                                }
                                result = true;
                            }
                        }
                        upliftVanSlidingCurve.RightPoint = new GeometryPoint(currentRightPoint.X, 0, currentRightPoint.Z);
                    }
                }
            }
            return result;
        }

        /// <summary>
        /// Determine the entry and exit points
        /// it is possible that there are more than 2 points
        /// Normaly this should be an even number
        /// The part chosen is with the greates didfference in Z level
        /// between entrance and exit point
        /// if the number of intersections point is odd 
        /// a valid intersection is found bij checking if the cirkel between two points lies insoil
        /// </summary>
        /// <returns></returns>
        private bool DetermineSurfacePoints(SlidingCircle slidingCircle)
        {
            // first check if intersection is not before or after begin or end surface;
            var intersectionpoints = LineHelper.ExtendedSurfaceIntersectionPointsWithCircle(slidingCircle.Center.X,
                                                                                            slidingCircle.Center.Y,
                                                                                            slidingCircle.Radius, surfaceLine.Points);
            intersectionpoints.Sort();
            if (intersectionpoints.Count == 2)
            {
                slidingCircle.LeftPoint = new GeometryPoint(intersectionpoints[0].X, 0, intersectionpoints[0].Z);
                slidingCircle.RightPoint = new GeometryPoint(intersectionpoints[1].X, 0, intersectionpoints[1].Z);
                return true;
            }

            if (intersectionpoints.Count > 1)
            {
                return HandleIntersectionsMoreThanTwo(slidingCircle, intersectionpoints);
            }
            return false;
        }

        private void DivideLargeSlicesInSmaller(SlidingCurve slidingCurve)
        {
            if (!double.IsNaN(striveWidth) && Math.Abs(striveWidth) > 0.01)
            {
                int numberOfSlices = slidingCurve.Slices.Count;

                for (int i = 0; i <= numberOfSlices - 2; i++)
                {
                    double dx = slidingCurve.Slices[i + 1].TopLeft.X - slidingCurve.Slices[i].TopLeft.X;
                    if (dx > striveWidth)
                    {
                        int numbers = Convert.ToInt32(Math.Floor(dx/striveWidth));
                        double size = dx/(numbers + 1);
                        for (int j = 1; j <= numbers; j++)
                        {
                            NextSlice(slidingCurve, slidingCurve.Slices[i].TopLeft.X + j*size);
                        }
                    }
                }
            }
        }

        private void addCentrePoints(SlidingDualCircle slidingCircle)
        {
            NextSlice(slidingCircle, slidingCircle.ActiveCircle.X);
            NextSlice(slidingCircle, slidingCircle.PassiveCircle.X);
        }

        private void addSurfacePoints(SlidingDualCircle slidingCircle)
        {
            NextSlice(slidingCircle, slidingCircle.LeftPoint.X);
            NextSlice(slidingCircle, slidingCircle.RightPoint.X);
        }

        private void addPoinstAboveVanSlipPlane(SlidingDualCircle slidingCircle)
        {
            // Deze punten zijn bij het preprocessen gemaakt
            double activeRadius = slidingCircle.ActiveCircle.Y - slidingCircle.TangentLine;
            double passiveRadius = slidingCircle.PassiveCircle.Y - slidingCircle.TangentLine;

            for (int i = 1; i < changingPoints.Count; i++)
            {
                // Only relevant geomery points
                if ((changingPoints[i - 1].X > slidingCircle.LeftPoint.X) &&
                    (changingPoints[i - 1].X < slidingCircle.RightPoint.X))
                {
                    if (
                        (Distance(changingPoints[i - 1].X, changingPoints[i - 1].Z, slidingCircle.ActiveCircle.X,
                                  slidingCircle.ActiveCircle.Y) <= activeRadius + Constants.CGeoAccu))
                    {
                        NextSlice(slidingCircle, changingPoints[i - 1].X);
                    }
                    else
                    {
                        if (
                            (Distance(changingPoints[i - 1].X, changingPoints[i - 1].Z,
                                      slidingCircle.PassiveCircle.X, slidingCircle.PassiveCircle.Y) <=
                             passiveRadius + Constants.CGeoAccu))
                        {
                            NextSlice(slidingCircle, changingPoints[i - 1].X);
                        }
                        else
                        {
                            if (changingPoints[i - 1].Y > slidingCircle.TangentLine)
                            {
                                if ((changingPoints[i - 1].X >
                                     Math.Min(slidingCircle.PassiveCircle.X, slidingCircle.ActiveCircle.X)) &&
                                    (changingPoints[i - 1].X <
                                     Math.Max(slidingCircle.PassiveCircle.X, slidingCircle.ActiveCircle.X)))
                                {
                                    NextSlice(slidingCircle, changingPoints[i - 1].X);
                                }
                            }
                        }
                    }
                }
            }
        }

        private void DetermineCircleData(SlidingDualCircle slidingCircle,
                                         ref double AXMid, ref double AZMid, ref double ARad, bool AIsLeftRequested)
        {
            if (AIsLeftRequested)
            {
                if (slidingCircle.LeftCircleIsActive)
                {
                    AXMid = slidingCircle.ActiveCircle.X;
                    AZMid = slidingCircle.ActiveCircle.Y;
                    ARad = slidingCircle.ActiveCircle.Y - slidingCircle.TangentLine;
                }
                else
                {
                    AXMid = slidingCircle.PassiveCircle.X;
                    AZMid = slidingCircle.PassiveCircle.Y;
                    ARad = slidingCircle.PassiveCircle.Y - slidingCircle.TangentLine;
                }
            }
            else
            {
                if (!slidingCircle.LeftCircleIsActive)
                {
                    AXMid = slidingCircle.ActiveCircle.X;
                    AZMid = slidingCircle.ActiveCircle.Y;
                    ARad = slidingCircle.ActiveCircle.Y - slidingCircle.TangentLine;
                }
                else
                {
                    AXMid = slidingCircle.PassiveCircle.X;
                    AZMid = slidingCircle.PassiveCircle.Y;
                    ARad = slidingCircle.PassiveCircle.Y - slidingCircle.TangentLine;
                }
            }
        }

        private void DetermineUpliftIntersectionPointsWithLayers(SlidingDualCircle slidingCircle)
        {
            double xMidLeft = 0;
            double xMidRight = 0;
            double zMidLeft = 0;
            double zMidRight = 0;
            double radLeft = 0;
            double radRight = 0;
            double xs = 0;
            double zs = 0;

            bool isLeftRequested = true;
            DetermineCircleData(slidingCircle, ref xMidLeft, ref zMidLeft, ref radLeft, isLeftRequested);

            isLeftRequested = false;
            DetermineCircleData(slidingCircle, ref xMidRight, ref zMidRight, ref radRight, isLeftRequested);

            foreach (var surface in profile2D.Surfaces)
            {
                bool finish = (surface.GeometrySurface.OuterLoop.GetMinX() > xMidRight + radRight);

                if (!finish)
                {
                    finish = surface.GeometrySurface.OuterLoop.GetMaxX() < xMidLeft - radLeft;
                }
                if (!finish)
                {
                    finish = surface.GeometrySurface.OuterLoop.GetMaxZ() < slidingCircle.TangentLine;
                }
                if (!finish)
                {
                    // Intersection is possible
                    for (int j = 0; j < surface.GeometrySurface.OuterLoop.Points.Count; j++)
                    {
                        int index = j == 0 ? surface.GeometrySurface.OuterLoop.Points.Count - 1 : j - 1;

                        double x1 = surface.GeometrySurface.OuterLoop.Points[index].X;
                        double z1 = surface.GeometrySurface.OuterLoop.Points[index].Z;
                        double x2 = surface.GeometrySurface.OuterLoop.Points[j].X;
                        double z2 = surface.GeometrySurface.OuterLoop.Points[j].Z;

                        int ns = 0;
                        double xs1 = 0;
                        double zs1 = 0;
                        double xs2 = 0;
                        double zs2 = 0;
                        
                        // See if there is intersection with left circel
                        MStabDatafunctions.Intersect_Circle_line(xMidLeft, zMidLeft, radLeft,
                                                                 x1, x2, z1, z2, ref ns, ref xs1, ref zs1,
                                                                 ref xs2, ref zs2);

                        if ((ns > 0))
                        {
                            if ((xs1 < xMidLeft))
                            {
                                NextSlice(slidingCircle, xs1);
                            }
                        }
                        if ((ns > 1))
                        {
                            if ((xs2 < xMidLeft))
                            {
                                NextSlice(slidingCircle, xs2);
                            }
                        }
                        // see if there is an intersection with the straight line
                        if (MStabDatafunctions.IntersectLines(x1, z1, x2, z2, xMidLeft,
                                                              slidingCircle.TangentLine, xMidRight,
                                                              slidingCircle.TangentLine, ref xs, ref zs))
                        {
                            NextSlice(slidingCircle, xs);
                        }
                        // See if there is intersection with the right circel
                        MStabDatafunctions.Intersect_Circle_line(xMidRight, zMidRight, radRight,
                                                                 x1, x2, z1, z2, ref ns, ref xs1, ref zs1,
                                                                 ref xs2, ref zs2);
                        if ((ns > 0))
                        {
                            if ((xs1 > xMidRight))
                            {
                                NextSlice(slidingCircle, xs1);
                            }
                        }
                        if ((ns > 1))
                        {
                            if ((xs2 > xMidRight))
                            {
                                NextSlice(slidingCircle, xs2);
                            }
                        }
                    }
                }
            }
        }

        private void DetermineUpliftIntersectionPointsWithWaternetlines(SlidingDualCircle slidingCircle)
        {
            double LXMidLeft = 0;
            double LXMidRight = 0;
            double LZMidLeft = 0;
            double LZMidRight = 0;
            double LRadLeft = 0;
            double LRadRight = 0;
            bool LIsLeftRequested = false;

            LIsLeftRequested = true;
            DetermineCircleData(slidingCircle, ref LXMidLeft, ref LZMidLeft, ref LRadLeft, LIsLeftRequested);
            LIsLeftRequested = false;
            DetermineCircleData(slidingCircle, ref LXMidRight, ref LZMidRight, ref LRadRight, LIsLeftRequested);

            foreach (var waternetLine in waternet.WaternetLineList)
            {
                AddSliceAtIntersection(slidingCircle, waternetLine, LXMidRight, LRadRight, LXMidLeft, LRadLeft, LZMidLeft, LZMidRight);
            }
            if (waternet.PhreaticLine != null)
            {
                AddSliceAtIntersection(slidingCircle, waternet.PhreaticLine, LXMidRight, LRadRight, LXMidLeft, LRadLeft, LZMidLeft, LZMidRight);
            }
        }

        private void AddSliceAtIntersection(SlidingDualCircle slidingCircle, GeometryPointString waternetLine, double LXMidRight,
                                            double LRadRight, double LXMidLeft, double LRadLeft, double LZMidLeft, double LZMidRight)

        {
            bool LFinish;
            int j;
            double Lx1;
            double LZ1;
            double LX2;
            double LZ2;
            double Lxs1 = 0;
            double LZs1 = 0;
            double LXs2 = 0;
            double LZs2 = 0;
            double LXS = 0;
            double LZS = 0;
            int LNs = 0;

            LFinish = (waternetLine.GetMinX() > LXMidRight + LRadRight);
            if (!LFinish)
            {
                LFinish = (waternetLine.GetMaxX() < LXMidLeft - LRadLeft);
            }
            if (!LFinish)
            {
                double MaxZ = waternetLine.GetMaxZ();
                LFinish = (MaxZ < slidingCircle.TangentLine);
            }
            if (!LFinish)
            {
                for (j = 1; j < waternetLine.Points.Count; j++)
                {
                    Lx1 = waternetLine.Points[j - 1].X;
                    LZ1 = waternetLine.Points[j - 1].Z;
                    LX2 = waternetLine.Points[j].X;
                    LZ2 = waternetLine.Points[j].Z;
                    // See if there is intersection with left circel
                    MStabDatafunctions.Intersect_Circle_line(LXMidLeft, LZMidLeft, LRadLeft,
                                                             Lx1, LX2, LZ1, LZ2, ref LNs, ref Lxs1, ref LZs1,
                                                             ref LXs2, ref LZs2);
                    if ((LNs > 0))
                    {
                        if ((Lxs1 < LXMidLeft))
                        {
                            NextSlice(slidingCircle, Lxs1);
                        }
                    }
                    if ((LNs > 1))
                    {
                        if ((LXs2 < LXMidLeft))
                        {
                            NextSlice(slidingCircle, LXs2);
                        }
                    }
                    // see if there is an intersection with the straight line
                    if (MStabDatafunctions.IntersectLines(Lx1, LZ1, LX2, LZ2, LXMidLeft,
                                                          slidingCircle.TangentLine, LXMidRight,
                                                          slidingCircle.TangentLine, ref LXS, ref LZS))
                    {
                        NextSlice(slidingCircle, LXS);
                    }
                    // See if there is intersection with the right circel
                    MStabDatafunctions.Intersect_Circle_line(LXMidRight, LZMidRight, LRadRight,
                                                             Lx1, LX2, LZ1, LZ2, ref LNs, ref Lxs1, ref LZs1,
                                                             ref LXs2, ref LZs2);
                    if ((LNs > 0))
                    {
                        if ((Lxs1 > LXMidRight))
                        {
                            NextSlice(slidingCircle, Lxs1);
                        }
                    }
                    if ((LNs > 1))
                    {
                        if ((LXs2 > LXMidRight))
                        {
                            NextSlice(slidingCircle, LXs2);
                        }
                    }
                }
            }
        }

        private void DetermineUpliftIntersectionPointsWithHeadlines(SlidingDualCircle slidingCircle)
        {
            int j = 0;
            bool LFinish = false;
            double Lx1 = 0;
            double LX2 = 0;
            double LZ1 = 0;
            double LZ2 = 0;
            int LNs = 0;
            double Lxs1 = 0;
            double LZs1 = 0;
            double LXs2 = 0;
            double LZs2 = 0;
            double LXMidLeft = 0;
            double LXMidRight = 0;
            double LZMidLeft = 0;
            double LZMidRight = 0;
            double LRadLeft = 0;
            double LRadRight = 0;
            bool LIsLeftRequested = false;
            double LXS = 0;
            double LZS = 0;

            LIsLeftRequested = true;
            DetermineCircleData(slidingCircle, ref LXMidLeft, ref LZMidLeft, ref LRadLeft, LIsLeftRequested);
            LIsLeftRequested = false;
            DetermineCircleData(slidingCircle, ref LXMidRight, ref LZMidRight, ref LRadRight, LIsLeftRequested);

            foreach (var headline  in waternet.HeadLineList)
            {
                LFinish = (headline.GetMinX() > LXMidRight + LRadRight);
                if (!LFinish)
                {
                    LFinish = (headline.GetMaxX() < LXMidLeft - LRadLeft);
                }
                if (!LFinish)
                {
                    LFinish = (headline.GetMaxZ() < slidingCircle.TangentLine);
                }
                if (!LFinish)
                {
                    for (j = 1; j < headline.Points.Count; j++)
                    {
                        Lx1 = headline.Points[j - 1].X;
                        LZ1 = headline.Points[j - 1].Z;
                        LX2 = headline.Points[j].X;
                        LZ2 = headline.Points[j].Z;
                        // See if there is intersection with left circel
                        MStabDatafunctions.Intersect_Circle_line(LXMidLeft, LZMidLeft, LRadLeft,
                                                                 Lx1, LX2, LZ1, LZ2, ref LNs, ref Lxs1, ref LZs1,
                                                                 ref LXs2, ref LZs2);
                        if ((LNs > 0))
                        {
                            if ((Lxs1 < LXMidLeft))
                            {
                                NextSlice(slidingCircle, Lxs1);
                            }
                        }
                        if ((LNs > 1))
                        {
                            if ((LXs2 < LXMidLeft))
                            {
                                NextSlice(slidingCircle, LXs2);
                            }
                        }
                        // see if there is an intersection with the straight line
                        if (MStabDatafunctions.IntersectLines(Lx1, LZ1, LX2, LZ2, LXMidLeft,
                                                              slidingCircle.TangentLine, LXMidRight,
                                                              slidingCircle.TangentLine, ref LXS, ref LZS))
                        {
                            NextSlice(slidingCircle, LXS);
                        }
                        // See if there is intersection with the right circel
                        MStabDatafunctions.Intersect_Circle_line(LXMidRight, LZMidRight, LRadRight,
                                                                 Lx1, LX2, LZ1, LZ2, ref LNs, ref Lxs1, ref LZs1,
                                                                 ref LXs2, ref LZs2);
                        if ((LNs > 0))
                        {
                            if ((Lxs1 > LXMidRight))
                            {
                                NextSlice(slidingCircle, Lxs1);
                            }
                        }
                        if ((LNs > 1))
                        {
                            if ((LXs2 > LXMidRight))
                            {
                                NextSlice(slidingCircle, LXs2);
                            }
                        }
                    }
                }
            } //========================================================================
        }
    }
}