using System;
using System.Collections.Generic;
using System.Linq;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Mathematics;
using Deltares.Standard.EventPublisher;
using Deltares.Standard.Language;
using Deltares.Standard.Logging;

namespace Deltares.Stability
{
    public class RTOCalculationgrid
    {
        private const bool isGridSearch = false;

        private GeometryPoint dikeShoulderLeft;
        private GeometryPoint dikeShoulderRight;
        private GeometryPoint ditchTop;

        private SoilProfile2D profile2D; // Aggregation
        private SurfaceLine2 surfaceLine; // Aggregation, not, is composition (see contracts)
        private Waternet waternet;

        public RTOCalculationgrid() { }

        /// <summary>
        /// Initializes a new instance of the <see cref="RTOCalculationgrid"/> class.
        /// </summary>
        /// <param name="stabilityModel">The stability model.</param>
        /// <Preconditions><paramref name="stabilityModel"/> cannot be null, 
        /// and its <see cref="StabilityModel.SoilProfile"/>, <see cref="StabilityModel.GeometryData"/>,
        /// <see cref="StabilityModel.SlipCircle"/> cannot be null.</Preconditions>
        /// <Postconditions><see cref="SurfaceLine"/> is not null.</Postconditions>
        public RTOCalculationgrid(StabilityModel stabilityModel)
        {
            profile2D = stabilityModel.SoilProfile;
            surfaceLine = stabilityModel.SurfaceLine2;
            waternet = stabilityModel.GeotechnicsData.CurrentWaternet;

            // remove the Nan
            RemoveNaNFromSurfaceLinePoints(surfaceLine);

            InwardStability = stabilityModel.SlipCircle.Orientation == GridOrientation.Inwards;
        }

        public bool InwardStability { get; set; }

        public SurfaceLine2 SurfaceLine // Aggregation, not at all: composition -> see contracts
        {
            get
            {
                return surfaceLine;
            }
            set
            {
                surfaceLine = value;
            }
        }

        /// <summary>
        /// Sets the soil profile2 d.
        /// </summary>
        /// <param name="profile">The profile.</param>
        /// <Precondtions><paramref name="profile"/> cannot be null.</Precondtions>
        public void SetSoilProfile2D(SoilProfile2D profile)
        {
            profile2D = profile;
        }

        /// <summary>
        /// Determines outward space to be filled with grids
        /// </summary>
        /// <param name="slipCircle"></param>
        /// <preconditions>
        /// <see cref="SurfaceLine"/> not null, its <see cref="BaseSurfaceLine.CharacteristicPoints"/> not null,
        /// and should contain characteristic points <see cref="CharacteristicPointType.DikeTopAtRiver"/> and 
        /// <see cref="CharacteristicPointType.DikeToeAtRiver"/>.
        /// </preconditions>
        public void DetermineOutwardBishopGrid(SlipCircle slipCircle)
        {
            GeometryPoint dikeTopRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeToeRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
            GeometryPoint dikeTopPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            double xMiddleDike = (dikeTopPolder.X + dikeTopRiver.X) / 2;

            slipCircle.SlipCircleGrid.GridXRight = dikeTopRiver.X;
            slipCircle.SlipCircleGrid.GridZBottom = dikeTopRiver.Z + 1.0;
            double extraheight = 0;
            if (IsOutsideShoulderPresent())
            {
                if (isGridSearch)
                {
                    slipCircle.SlipCircleGrid.GridZBottom = dikeShoulderLeft.Z;
                    extraheight = dikeTopRiver.Z - dikeShoulderLeft.Z;
                }
                else
                {
                    slipCircle.SlipCircleGrid.GridZBottom = dikeShoulderLeft.Z;
                }
            }

            slipCircle.SlipCircleGrid.GridXLeft = dikeToeRiver.X;

            double height = dikeTopRiver.Z - GetTopPleistoceen(xMiddleDike) + extraheight;
            slipCircle.SlipCircleGrid.GridZTop = slipCircle.SlipCircleGrid.GridZBottom + height;

            slipCircle.SlipCircleTangentLine.TangentLineZTop = dikeToeRiver.Z + 2.0;
            slipCircle.SlipCircleTangentLine.TangentLineZBottom = Math.Min(GetTopPleistoceen(xMiddleDike), dikeToeRiver.Z);
        }

        /// <summary>
        /// Create the grid boundaries for a grid calculation with search minimum
        /// </summary>
        /// <param name="slipCircle"> </param>
        /// <preconditions>
        /// <see cref="SurfaceLine"/> cannot be null, its <see cref="BaseSurfaceLine.CharacteristicPoints"/> cannot be null
        /// and <see cref="GeometryPointString.Points"/> should have points.
        /// </preconditions>
        public void UpdateBishopGrid(SlipCircle slipCircle)
        {
            if (slipCircle == null)
            {
                throw new ArgumentNullException("slipCircle");
            }
            if (SurfaceLine == null)
            {
                throw new InvalidOperationException("Grid can not be generated / updated when surface line is not present.");
            }
            
            if (InwardStability)
            {
                DetermineInwardBishopGrid(slipCircle);
            }
            else
            // outward stability = river side
            {
                DetermineOutwardBishopGrid(slipCircle);
            }

            // determine the intervalls
            slipCircle.SlipCircleGrid.XIntervalNumber = 0;
            if (Math.Abs(slipCircle.SlipCircleGrid.GridXLeft - slipCircle.SlipCircleGrid.GridXRight) >
                Constants.CGeoAccu)
            {
                slipCircle.SlipCircleGrid.XIntervalNumber =
                    (int)Math.Max(Math.Round(Math.Abs(slipCircle.SlipCircleGrid.GridXRight - slipCircle.SlipCircleGrid.GridXLeft)), 2);
            }
            slipCircle.SlipCircleGrid.ZIntervalNumber = 0;
            if (Math.Abs(slipCircle.SlipCircleGrid.GridZTop - slipCircle.SlipCircleGrid.GridZBottom) >
                Constants.CGeoAccu)
            {
                slipCircle.SlipCircleGrid.ZIntervalNumber =
                    (int)Math.Max(Math.Round(Math.Abs(slipCircle.SlipCircleGrid.GridZTop - slipCircle.SlipCircleGrid.GridZBottom)), 2);
            }
            slipCircle.SlipCircleTangentLine.TangentLineNumber = 0;
            if (
                Math.Abs(slipCircle.SlipCircleTangentLine.TangentLineZTop -
                         slipCircle.SlipCircleTangentLine.TangentLineZBottom) > Constants.CGeoAccu)
            {
                slipCircle.SlipCircleTangentLine.TangentLineNumber =
                    (int)
                    Math.Max(
                        Math.Round(
                            Math.Abs(slipCircle.SlipCircleTangentLine.TangentLineZTop -
                                     slipCircle.SlipCircleTangentLine.TangentLineZBottom)), 3);
            }
        }

        public void CreateZone(StabilityModel model)
        {
            model.SlipPlaneConstraints.XEntryMin = double.NaN;
            model.SlipPlaneConstraints.XEntryMax = double.NaN;

            if (surfaceLine.IsDefined(CharacteristicPointType.DikeToeAtRiver))
            {
                model.SlipPlaneConstraints.XEntryMin = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X;
            }

            if (surfaceLine.IsDefined(CharacteristicPointType.DikeTopAtPolder)
                && surfaceLine.IsDefined(CharacteristicPointType.DikeToeAtPolder))
            {
                GeometryPoint top = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
                GeometryPoint toe = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);

                double zEnd = 0.5 * (top.Z + toe.Z);
                if (surfaceLine.IsDefined(CharacteristicPointType.ShoulderBaseInside))
                {
                    GeometryPoint shoulder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);

                    if (zEnd < shoulder.Z)
                    {
                        model.SlipPlaneConstraints.XEntryMax = shoulder.X;
                    }
                    else
                    {
                        model.SlipPlaneConstraints.XEntryMax = FindXAtSurface(top.X, toe.X, zEnd);
                    }
                }
                else
                {
                    model.SlipPlaneConstraints.XEntryMax = FindXAtSurface(top.X, toe.X, zEnd);
                }
            }
        }

        /// <summary>
        /// Creates a calculation grid if autogeneration is requested.
        /// </summary>
        /// <param name="model">Model with data to use for generating the grid.</param>
        /// <exception cref="ArgumentException">if <paramref name="model"/> has no soil layers and grid should be generated.</exception>
        public void CreateGrid(StabilityModel model)
        {
            if (!model.SoilProfile.Surfaces.Any())
            {
                return;
            }

            bool autoGenerate = false;
            switch (model.ModelOption)
            {
                case ModelOptions.Bishop:
                case ModelOptions.UpliftVan:
                    autoGenerate = model.SlipCircle.Auto;
                    break;
                case ModelOptions.Spencer:
                    autoGenerate = model.AutoGenerateGeneticSpencer;
                    break;
            }

            if (autoGenerate)
            {
                waternet = model.GeotechnicsData.CurrentWaternet;

                if (StabilityModelHasLayers(model))
                {
                    if (model.SearchAlgorithm == SearchAlgorithm.Genetic ||
                        model.SearchAlgorithm == SearchAlgorithm.LevenbergMarquardt ||
                        model.SearchAlgorithm == SearchAlgorithm.GeneticAndLevenbergMarquardt)
                    {
                        switch (model.ModelOption)
                        {
                            case ModelOptions.Spencer:
                                UpdateSpencerGA(model);
                                break;
                        }
                    }
                    else if (model.SearchAlgorithm == SearchAlgorithm.Grid)
                    {
                        switch (model.ModelOption)
                        {
                            case ModelOptions.Bishop:
                                UpdateBishopGrid(model.SlipCircle);
                                break;
                            case ModelOptions.UpliftVan:
                                UpdateUpliftGrid(model.SlipPlaneUpliftVan);
                                break;
                            case ModelOptions.Spencer:
                                UpdateSpencerGA(model);
                                break;
                        }
                    }
                }
            }
        }

        public void GenerateSpencerSlipPlanes(StabilityModel model)
        {
            if (model.IsHighSlipPlanes)
            {
                CreateInwardSpencerPlanesHigh(model);
            }
            else
            {
                CreateInwardSpencerPlanesLow(model);
            }
        }

        public void ReplaceSlipPlanesObject(object target)
        {
            if (target is StabilityModel)
            {
                ReplaceSlipPlanes(target as StabilityModel);
            }
        }

        public void ReplaceSlipPlanes(StabilityModel stabilityModel)
        {
            DataEventPublisher.BeforeChange(this, "SlipPlanes");

            ResetSlipPlanes(stabilityModel);

            if (stabilityModel.SurfaceLine2 != null)
            {
                DataEventPublisher.InvokeWithoutPublishingEvents(() =>
                    {
                        try
                        {
                            var rtoCalculationGrid = new RTOCalculationgrid(stabilityModel);
                            rtoCalculationGrid.CreateGrid(stabilityModel);

                            foreach (var slipPlane in stabilityModel.SlipPlanes)
                            {
                                slipPlane.Name = LocalizationManager.GetTranslatedText(typeof(StabilityDikeLocationInfo), "SlipPlane");
                            }
                        }
                        catch (Exception)
                        {
                            ResetSlipPlanes(stabilityModel);

                            if (stabilityModel.SurfaceLine2.CharacteristicPoints.Any() && stabilityModel.SoilProfile.Surfaces.Any()) // Only add log message when no log message is already generated via validation
                            {
                                var logMessage = new LogMessage(LogMessageType.Warning, this, string.Format(LocalizationManager.GetTranslatedText(typeof(StabilityDikeLocationInfo), "CannotCreateSlipPlanes")));

                                if (LogManager.Messages.All(m => m.Message != logMessage.Message)) // Only add log message once
                                {
                                    LogManager.Add(logMessage);
                                }
                            }
                        }
                    });
            }

            DataEventPublisher.AfterChange(this, "SlipPlanes");
        }

        /// <summary>
        /// Find level from pleistocene at X.
        /// If no aquitard is found, it is assumed that the bottom layer is the pleistocene level.
        /// This is done by creating a 1D geometry at point X. Take the bottom level as first guess.
        /// </summary>
        /// <preconditions><see cref="profile2D"/> or <see cref="surfaceLine"/> are not null, and <see cref="surfaceLine"/>
        /// has the following characteristic points: <see cref="CharacteristicPointType.DikeTopAtRiver"/> and
        /// <see cref="CharacteristicPointType.DikeTopAtPolder"/>. 1D profile corresponding to <see cref="profile2D"/>
        /// for the given characteristic points has layers.</preconditions>
        private double GetTopPleistoceen(double x)
        {
            SoilProfile1D helpProfile = profile2D.GetSoilProfile1D(x);

            if (helpProfile.Layers.Count == 0)
            {
                return double.NaN;
            }

            double result = helpProfile.Layers[helpProfile.LayerCount - 1].TopLevel;
            for (int i = helpProfile.LayerCount - 2; i > -1; i--)
            {
                if (!helpProfile.Layers[i].IsAquifer)
                {
                    result = helpProfile.Layers[i + 1].TopLevel;
                    break;
                }
            }
            return result;
        }

        /// <summary>
        /// Fins out if there is a shoulder is knwn at the polder side
        /// </summary>
        /// <returns></returns>
        private bool IsInsideShoulderPresent()
        {
            bool result = false;
            dikeShoulderLeft = surfaceLine.IsDefined(CharacteristicPointType.ShoulderBaseInside)
                                   ? surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside)
                                   : null;
            dikeShoulderRight = surfaceLine.IsDefined(CharacteristicPointType.ShoulderTopInside)
                                    ? surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderTopInside)
                                    : null;
            if ((dikeShoulderLeft != null) && (dikeShoulderRight != null))
            {
                result = (dikeShoulderLeft.X < dikeShoulderRight.X);
            }

            return result;
        }

        /// <summary>
        /// find out if a ditch is given (only at polder side)
        /// </summary>
        /// <returns></returns>
        private bool IsDitchPresent()
        {
            bool result = false;
            ditchTop = surfaceLine.IsDefined(CharacteristicPointType.DitchPolderSide)
                           ? surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide)
                           : null;
            GeometryPoint ditchBot = surfaceLine.IsDefined(CharacteristicPointType.BottomDitchPolderSide)
                                         ? surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchPolderSide)
                                         : null;
            if ((ditchTop != null) && (ditchBot != null))
            {
                if ((!double.IsNaN(ditchTop.X)) && (!double.IsNaN(ditchBot.X)))
                {
                    result = true;
                }
            }
            return result;
        }

        /// <summary>
        /// find out if a ditch is given (only at polder side)
        /// </summary>
        /// <returns></returns>
        private bool IsChannelPresent()
        {
            return surfaceLine.IsDefined(CharacteristicPointType.InsertRiverChannel);
        }

        /// <summary>
        /// Check if a shoulder is fiven at de river side
        /// </summary>
        /// <returns></returns>
        private bool IsOutsideShoulderPresent()
        {
            bool result = false;
            dikeShoulderLeft = surfaceLine.IsDefined(CharacteristicPointType.ShoulderTopOutside)
                                   ? surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderTopOutside)
                                   : null;
            dikeShoulderRight = surfaceLine.IsDefined(CharacteristicPointType.ShoulderBaseOutside)
                                    ? surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseOutside)
                                    : null;
            if ((dikeShoulderLeft != null) && (dikeShoulderRight != null))
            {
                result = (dikeShoulderLeft.X < dikeShoulderRight.X);
            }
            return result;
        }

        /// <summary>
        /// Determines inward space to be filled with grids
        /// </summary>
        /// <param name="slipCircle"></param>
        /// <preconditions>
        /// <see cref="SurfaceLine"/> not null, its <see cref="BaseSurfaceLine.CharacteristicPoints"/> not null,
        /// and should contain characteristic points <see cref="CharacteristicPointType.DikeTopAtPolder"/> and 
        /// <see cref="CharacteristicPointType.DikeToeAtPolder"/>.
        /// </preconditions>
        private void DetermineInwardBishopGrid(SlipCircle slipCircle)
        {
            GeometryPoint dikeTopPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint dikeToe = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint dikeTopRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            double xMiddleDike = (dikeTopPolder.X + dikeTopRiver.X) / 2;

            slipCircle.SlipCircleGrid.GridXLeft = dikeTopPolder.X;
            slipCircle.SlipCircleGrid.GridZBottom = dikeTopPolder.Z + 1.0;
            double extraheight = 0;
            if (IsInsideShoulderPresent())
            {
                if (isGridSearch)
                {
                    slipCircle.SlipCircleGrid.GridZBottom = dikeShoulderRight.Z;
                    extraheight = dikeTopPolder.Z - dikeShoulderRight.Z;
                }
                else
                {
                    slipCircle.SlipCircleGrid.GridZBottom = dikeShoulderRight.Z;
                }
            }

            double height = dikeTopPolder.Z - GetTopPleistoceen(xMiddleDike) + extraheight;
            slipCircle.SlipCircleGrid.GridZTop = slipCircle.SlipCircleGrid.GridZBottom + height;
            slipCircle.SlipCircleGrid.GridXRight = dikeToe.X;
            slipCircle.SlipCircleTangentLine.TangentLineZTop = dikeToe.Z + 2.0;
            slipCircle.SlipCircleTangentLine.TangentLineZBottom = Math.Min(GetTopPleistoceen(xMiddleDike), dikeToe.Z);
        }

        /// <summary>
        /// Determine space to be filled with grid (only polderside)
        /// </summary>
        /// <param name="slipPlaneUpliftVan"></param>
        private void UpdateUpliftGrid(SlipPlaneUpliftVan slipPlaneUpliftVan)
        {
            InwardStability = true;
            // ACTIVE GRID is the same as Bishop
            var tempCircle = new SlipCircle();
            UpdateBishopGrid(tempCircle);
            slipPlaneUpliftVan.SlipPlaneLeftGrid = tempCircle.SlipCircleGrid;
            slipPlaneUpliftVan.SlipPlaneLeftGrid.XIntervalNumber = (int)Math.Max(
                Math.Round(Math.Abs((slipPlaneUpliftVan.SlipPlaneLeftGrid.GridXRight -
                                     slipPlaneUpliftVan.SlipPlaneLeftGrid.GridXLeft) / 1.5)), 2);
            slipPlaneUpliftVan.SlipPlaneLeftGrid.ZIntervalNumber = (int)Math.Max(
                Math.Round(Math.Abs((slipPlaneUpliftVan.SlipPlaneLeftGrid.GridZTop -
                                     slipPlaneUpliftVan.SlipPlaneLeftGrid.GridZBottom) / 1.5)), 2);

            // PASSIVE GRID 
            GeometryPoint dikeToe = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint rightBoundary = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside);
            GeometryPoint ditchDikeSide = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);

            slipPlaneUpliftVan.SlipPlaneRightGrid.GridXLeft = dikeToe.X;

            if (IsDitchPresent())
            {
                slipPlaneUpliftVan.SlipPlaneRightGrid.GridXRight = ditchDikeSide.X;
                slipPlaneUpliftVan.SlipPlaneRightGrid.GridZBottom = ditchDikeSide.Z;
            }
            else
            {
                double xCoordAtWhichUpliftEnds = SearchXcoordinateOfEndOfUpliftAtPolderSide();
                // If no uplift occurs, Xright of passive grid is the dike toe + 8 m
                if (Math.Abs(xCoordAtWhichUpliftEnds - rightBoundary.X) < Constants.CGeoAccu)
                {
                    slipPlaneUpliftVan.SlipPlaneRightGrid.GridXRight = Math.Min(dikeToe.X + 8.0, rightBoundary.X);
                }
                else
                {
                    slipPlaneUpliftVan.SlipPlaneRightGrid.GridXRight = xCoordAtWhichUpliftEnds;
                }
                slipPlaneUpliftVan.SlipPlaneRightGrid.GridZBottom = dikeToe.Z;
            }

            slipPlaneUpliftVan.SlipPlaneRightGrid.GridZTop = slipPlaneUpliftVan.SlipPlaneRightGrid.GridZBottom + 4.0;

            slipPlaneUpliftVan.SlipPlaneRightGrid.XIntervalNumber = 0;
            if (Math.Abs(slipPlaneUpliftVan.SlipPlaneRightGrid.GridXLeft - slipPlaneUpliftVan.SlipPlaneRightGrid.GridXRight) > Constants.CGeoAccu)
            {
                slipPlaneUpliftVan.SlipPlaneRightGrid.XIntervalNumber = (int)Math.Max(Math.Round((slipPlaneUpliftVan.SlipPlaneRightGrid.GridXRight -
                                                                                                   slipPlaneUpliftVan.SlipPlaneRightGrid.GridXLeft) / 2.0), 2);
            }
            slipPlaneUpliftVan.SlipPlaneRightGrid.ZIntervalNumber = 0;
            if (Math.Abs(slipPlaneUpliftVan.SlipPlaneRightGrid.GridZTop - slipPlaneUpliftVan.SlipPlaneRightGrid.GridZBottom) > Constants.CGeoAccu)
            {
                slipPlaneUpliftVan.SlipPlaneRightGrid.ZIntervalNumber = 2;
            }

            // TANGENT LINES
            slipPlaneUpliftVan.SlipCircleTangentLine.AutomaticAtBoundaries = true;
            slipPlaneUpliftVan.SlipCircleTangentLine.MaxSpacingBetweenBoundaries = 1;

            slipPlaneUpliftVan.SlipPlaneTangentLine.UpdateBoundaries();
            slipPlaneUpliftVan.SlipCircleTangentLine.TangentLineZBottom = slipPlaneUpliftVan.SlipCircleTangentLine.BoundaryHeights.Max().Height;
            slipPlaneUpliftVan.SlipCircleTangentLine.TangentLineZTop = slipPlaneUpliftVan.SlipCircleTangentLine.BoundaryHeights.Min().Height;
        }

        private void RemoveNaNFromSurfaceLinePoints(SurfaceLine2 line)
        {
            if (line != null)
            {
                var pointsToRemove = line.Geometry.Points.Where(p => double.IsNaN(p.X) || double.IsNaN(p.Z)).ToList();
                foreach (var geometryPoint in pointsToRemove)
                {
                    line.Geometry.Points.Remove(geometryPoint);
                }
            }
        }

        private double FindXAtSurface(double startX, double endX, double level)
        {
            double result = Double.NaN;
            for (int i = 0; i < surfaceLine.Geometry.Count - 2; i++)
            {
                if ((surfaceLine.Geometry.Points[i].X >= startX) && (surfaceLine.Geometry.Points[i + 1].X <= endX))
                {
                    if ((Math.Max(surfaceLine.Geometry.Points[i].Z, surfaceLine.Geometry.Points[i + 1].Z) >= level) &&
                        (Math.Min(surfaceLine.Geometry.Points[i].Z, surfaceLine.Geometry.Points[i + 1].Z) <= level))
                    {
                        result = Routines2D.LinInpolY(surfaceLine.Geometry.Points[i].Z, surfaceLine.Geometry.Points[i].X,
                                                      surfaceLine.Geometry.Points[i + 1].Z, surfaceLine.Geometry.Points[i + 1].X, level);
                        break;
                    }
                }
            }
            return result;
        }

        private static bool StabilityModelHasLayers(StabilityModel model)
        {
            return model.SoilProfile.Surfaces.Any();
        }

        private void CreateInwardSpencerInnerSlipPlaneLow(GeometryPointString innerSlipPlane)
        {
            innerSlipPlane.Points.Clear();

            var entryPoint = new GeometryPoint();
            var exitPoint = new GeometryPoint();
            GeometryPoint dikeTopPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint dikeToe = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint ditchDikeSide = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);

            if (IsInsideShoulderPresent())
            {
                dikeShoulderLeft = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside);
                entryPoint.X = 0.5 * (dikeTopPolder.X + dikeShoulderLeft.X);
                entryPoint.Z = surfaceLine.Geometry.GetZAtX(entryPoint.X);
            }
            else
            {
                GeometryPoint diketoe = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
                entryPoint.X = 0.5 * (dikeTopPolder.X + diketoe.X);
                entryPoint.Z = surfaceLine.Geometry.GetZAtX(entryPoint.X);
            }
            innerSlipPlane.Points.Add(entryPoint);

            // Exit point (9in) = rightest point
            if (IsDitchPresent())
            {
                exitPoint.X = Math.Min(dikeToe.X + 5, ditchDikeSide.X);
                exitPoint.Z = surfaceLine.Geometry.GetZAtX(exitPoint.X);
            }
            else
            {
                exitPoint.X = dikeToe.X + 5;
                exitPoint.Z = surfaceLine.Geometry.GetZAtX(exitPoint.X);
            }

            // Height h is calculated at exit point (9in)
            double h = exitPoint.Z - GetTopPleistoceen(exitPoint.X);

            // First determine if point 6in is situated after the exit point. In such case, point 6in is the middle of the entry and exit point.
            bool point6inIsAtRightSideOfExitPoint = entryPoint.X + entryPoint.Z - (exitPoint.Z - h / 2) > exitPoint.X;
            var point6in = new GeometryPoint();
            if (point6inIsAtRightSideOfExitPoint)
            {
                point6in.X = (entryPoint.X + exitPoint.X) / 2;
            }
            else
            {
                point6in.X = entryPoint.X + entryPoint.Z - (exitPoint.Z - h / 2);
            }
            point6in.Z = exitPoint.Z - h / 2;

            // Then determine if left and right slopes (1:1) of the inner plane intersect above level Z_exit - h/2
            bool bothSlopesIntersect = point6in.X > exitPoint.X - h / 2;

            // Case 1: left and right slopes (1:1) of the inner plane do not intersect above level Z_exit - h/2
            if (!bothSlopesIntersect)
            {
                double deltaX = point6in.X - entryPoint.X;
                double deltaZ = entryPoint.Z - point6in.Z;
                // the 6 points on the left slope 1:1
                for (int i = 1; i < 6; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = entryPoint.X + i * deltaX / 5;
                    nextPoint.Z = entryPoint.Z - i * deltaZ / 5;
                    innerSlipPlane.Points.Add(nextPoint);
                }
                // point 8in
                var point8in = new GeometryPoint();
                if (h < Constants.CAlmostZero)
                {
                    point8in.X = exitPoint.X - (exitPoint.X - point6in.X) / 3;
                    point8in.Z = exitPoint.Z;
                }
                else
                {
                    point8in.X = exitPoint.X - h / 2;
                    point8in.Z = exitPoint.Z - h / 2;
                }
                // point 7in = middle of point 6in and 8 in
                var point7in = new GeometryPoint();
                point7in.X = (entryPoint.X + deltaX + point8in.X) / 2;
                point7in.Z = point8in.Z;
                innerSlipPlane.Points.Add(point7in);
                innerSlipPlane.Points.Add(point8in);
            }
            // Case 2: left and right slopes (1:1) of the hight inner plane do intersect
            else
            {
                // point 7in = intersection point of left and right slopes
                var point7in = new GeometryPoint();
                point7in.X = (entryPoint.Z + entryPoint.X - exitPoint.Z + exitPoint.X) / 2;
                point7in.Z = exitPoint.Z - (exitPoint.X - point7in.X);

                // point 8in = between 7in and 9in
                var point8in = new GeometryPoint();
                point8in.X = (exitPoint.X + point7in.X) / 2;
                point8in.Z = (exitPoint.Z + point7in.Z) / 2;

                // the 6 points (2in until 7in) on the left slope 1:1 are added
                double delta = point7in.X - entryPoint.X;
                for (int i = 1; i < 7; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = entryPoint.X + i * delta / 6;
                    nextPoint.Z = entryPoint.Z - i * delta / 6;
                    innerSlipPlane.Points.Add(nextPoint);
                }
                innerSlipPlane.Points.Add(point8in);
            }
            // exit point 9in
            innerSlipPlane.Points.Add(exitPoint);

            innerSlipPlane = EnsureThatPointsOfTheSlipPlaneAreNotAboveSurfaceLine(innerSlipPlane);

            DataEventPublisher.DataListModified(innerSlipPlane.Points);
        }

        private void CreateInwardSpencerInnerSlipPlaneHigh(GeometryPointString innerSlipPlane)
        {
            innerSlipPlane.Points.Clear();

            var exitPoint = new GeometryPoint();
            var entryPoint = new GeometryPoint();
            GeometryPoint dikeTopAtPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint dikeToeAtPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint dikeTopAtRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            double height = (dikeTopAtPolder.Z - dikeToeAtPolder.Z) / 3;

            // Exit point (9in) = rightest point
            exitPoint.X = dikeToeAtPolder.X;
            exitPoint.Z = dikeToeAtPolder.Z;

            // Entry point (1in) = leftest point 
            entryPoint.Z = exitPoint.Z + height;
            entryPoint.X = surfaceLine.Geometry.GetXAtZ(entryPoint.Z, dikeTopAtRiver.X);

            // Intersection point 
            var intersectionPoint = new GeometryPoint();
            intersectionPoint.X = entryPoint.X + (entryPoint.Z - exitPoint.Z);
            intersectionPoint.Z = Math.Min(surfaceLine.Geometry.GetZAtX(intersectionPoint.X), exitPoint.Z);

            innerSlipPlane.Points.Add(entryPoint);

            // Point 7in
            var point7in = new GeometryPoint();
            if (intersectionPoint.X > exitPoint.X)
            {
                point7in.X = (entryPoint.X + exitPoint.X) / 2;
                point7in.Z = Math.Min(surfaceLine.Geometry.GetZAtX(point7in.X), exitPoint.Z);
            }
            else
            {
                point7in = intersectionPoint;
            }

            // the seven points on the line 1:1
            for (int i = 1; i < 7; i++)
            {
                var nextPoint = new GeometryPoint();
                nextPoint.X = entryPoint.X + i / 6.0 * Math.Abs(point7in.X - entryPoint.X);
                nextPoint.Z = entryPoint.Z - i / 6.0 * Math.Abs(point7in.Z - entryPoint.Z);
                innerSlipPlane.Points.Add(nextPoint);
            }

            // Point 8in
            var point8in = new GeometryPoint();
            point8in.X = (point7in.X + exitPoint.X) / 2;
            point8in.Z = Math.Min(surfaceLine.Geometry.GetZAtX(point8in.X), exitPoint.Z);

            innerSlipPlane.Points.Add(point8in);
            innerSlipPlane.Points.Add(exitPoint);

            innerSlipPlane = EnsureThatPointsOfTheSlipPlaneAreNotAboveSurfaceLine(innerSlipPlane);

            DataEventPublisher.DataListModified(innerSlipPlane.Points);
        }

        /*
         
           start (point 1)   end (point 10)
                |            /                             |            /                        |            /
                |           /                              |           /                         |           /
                |          /                               |          /                          |          /
             1  |      2  /                                |         /                           |         /
                |        /                                 |        /                            |        /
                |       /                                  |       /                             |       /
                |      /                                   |      /                           -------------  bottom pleistoceen
                |     /                                    |     /                               |     /
                |    /                                     |    /                                |    /
                |   /                                      |   /                                 |   /
                |  /                                   ------------    bottom pleistoceen        |  /
                | /                                        | /                                   | /
                |/                                         |/                                    |/
             intersection

        */

        private void CreateInwardSpencerOuterSlipPlaneLow(GeometryPointString outerSlipPlane)
        {
            outerSlipPlane.Points.Clear();
            bool isDitchPresent = IsDitchPresent();

            GeometryPoint dikeTopRiver =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeToePolder =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint dikeTopPolder =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint ditchDikeSide =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
            GeometryPoint bottomDitchDikeSide =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchDikeSide);
            GeometryPoint bottomDitchPolderSide =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchPolderSide);
            GeometryPoint ditchPolderSide =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide);
            GeometryPoint surfaceLevelInside =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside);

            // Height h is calculated at exit point of the inner plane (Low) : 9in = min (N; M + 5m)
            var exitPointInnerPlane = new GeometryPoint();
            double d = 0;
            if (isDitchPresent)
            {
                exitPointInnerPlane.X = Math.Min(dikeToePolder.X + 5, ditchDikeSide.X);
                exitPointInnerPlane.Z = surfaceLine.Geometry.GetZAtX(exitPointInnerPlane.X);
                d = (ditchDikeSide.Z + ditchPolderSide.Z) / 2 - (bottomDitchDikeSide.Z + bottomDitchPolderSide.Z) / 2;
            }
            else
            {
                exitPointInnerPlane.X = dikeToePolder.X + 5;
                exitPointInnerPlane.Z = surfaceLine.Geometry.GetZAtX(exitPointInnerPlane.X);
            }
            double h = exitPointInnerPlane.Z - GetTopPleistoceen(exitPointInnerPlane.X);
            double bottom;
            if (h < 2) // Case 3
            {
                bottom = dikeToePolder.Z - 5.0;
            }
            else // Cases 1 and 2
            {
                bottom = GetTopPleistoceen(exitPointInnerPlane.X) - 1;
            }
            // Determine the exit point (9out) 
            var exitPoint = new GeometryPoint();
            if (h < 2) // case 3
            {
                exitPoint.X = Math.Min(dikeToePolder.X + 32.5, surfaceLevelInside.X - 0.01);
                exitPoint.Z = surfaceLine.Geometry.GetZAtX(exitPoint.X);
            }
            else // h > 2 m
            {
                if ((!isDitchPresent) || (isDitchPresent && (h / d > 5))) // Case 2 (shallow ditch or no ditch)
                {
                    exitPoint.X = Math.Min(dikeToePolder.X + 5.5 * (h + 1), surfaceLevelInside.X - 0.01);
                    exitPoint.Z = surfaceLine.Geometry.GetZAtX(exitPoint.X);
                }
                else // case 1: deep ditch (h / d < 5)
                {
                    exitPoint.X = bottomDitchPolderSide.X;
                    exitPoint.Z = bottomDitchPolderSide.Z;
                }
            }

            // Slope of the exit line
            double finalSlope;
            if (h < 2) // case 3
            {
                finalSlope = 2.5;
            }
            else
            {
                if (h / d < 5) // case 1
                {
                    finalSlope = 1;
                }
                else // case 2
                {
                    finalSlope = 2.5;
                }
            }

            // Determine the intersection of the final slope and the vertical part
            var intersection = new GeometryPoint();
            intersection.X = dikeTopRiver.X;
            intersection.Z = exitPoint.Z - (exitPoint.X - dikeTopRiver.X) / finalSlope;
            bool isDeepPleistocene = intersection.Z >= bottom;
            double vertLength;
            if (!isDeepPleistocene)
            {
                vertLength = dikeTopRiver.Z - bottom;
            }
            else
            {
                vertLength = dikeTopRiver.Z - intersection.Z;
            }
            var vertSegmentLength = vertLength / 4;

            // Entry point (1out) = intersection with line under 45 degrees starting at point 2out
            var fictivePoint1 = new GeometryPoint(); // = point 2out
            fictivePoint1.X = dikeTopRiver.X;
            fictivePoint1.Z = dikeTopRiver.Z - vertSegmentLength;

            var fictivePoint2 = new GeometryPoint();
            fictivePoint2.X = fictivePoint1.X;
            fictivePoint2.Z = fictivePoint1.Z;

            while (fictivePoint2.Z <= surfaceLine.Geometry.GetMaxZ())
            {
                fictivePoint2.X = fictivePoint2.X - vertSegmentLength;
                fictivePoint2.Z = fictivePoint2.Z + vertSegmentLength;
            }

            var obliqueLine = new Line(fictivePoint1, fictivePoint2);
            GeometryPoint entryPoint = DetermineIntersectionBetweenDikeSlopeAtRiverAndAnObliqueLine(surfaceLine, obliqueLine);
            outerSlipPlane.Points.Add(entryPoint);

            if (!isDeepPleistocene)
            {
                // The 3 points on the vertical line    
                for (int i = 1; i < 4; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = dikeTopRiver.X;
                    nextPoint.Z = dikeTopRiver.Z - i * vertSegmentLength;
                    outerSlipPlane.Points.Add(nextPoint);
                }

                // Point 8out
                var point8out = new GeometryPoint();
                point8out.Z = bottom;
                point8out.X = exitPoint.X - finalSlope * (exitPoint.Z - point8out.Z);

                // Point 6out
                var point6out = new GeometryPoint();
                var horLength = point8out.X - dikeTopRiver.X;
                var horSegmLength = horLength / 3;
                point6out.X = dikeTopRiver.X + horSegmLength;
                point6out.Z = point8out.Z;

                // Point 5out
                var point5out = new GeometryPoint();
                point5out.X = (point6out.X + dikeTopRiver.X) / 2;
                point5out.Z = (point6out.Z + dikeTopRiver.Z - 3.0 * vertSegmentLength) / 2;

                // Point 7out
                var point7out = new GeometryPoint();
                point7out.X = point6out.X + horSegmLength;
                point7out.Z = point8out.Z;

                outerSlipPlane.Points.Add(point5out);
                outerSlipPlane.Points.Add(point6out);
                outerSlipPlane.Points.Add(point7out);
                outerSlipPlane.Points.Add(point8out);
                outerSlipPlane.Points.Add(exitPoint);
            }
            // Other configuration (i.e. deep Pleistocene): Outer slip plane has no horizontal bottom part, 
            // but only a first vertical part and an oblique part
            else
            {
                // The 4 points on the vertical line (2out to 5out)
                for (int i = 1; i < 5; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = dikeTopRiver.X;
                    nextPoint.Z = dikeTopRiver.Z - i * (dikeTopRiver.Z - intersection.Z) / 4;
                    outerSlipPlane.Points.Add(nextPoint);
                }

                // The 4 points on the final slope (6out to 9out)
                for (int i = 1; i < 5; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = intersection.X + i * (exitPoint.X - intersection.X) / 4;
                    nextPoint.Z = intersection.Z + i * (exitPoint.Z - intersection.Z) / 4;
                    outerSlipPlane.Points.Add(nextPoint);
                }
            }

            outerSlipPlane = EnsureThatPointsOfTheSlipPlaneAreNotAboveSurfaceLine(outerSlipPlane);
            DataEventPublisher.DataListModified(outerSlipPlane.Points);
        }

        private void CreateInwardSpencerOuterSlipPlaneHigh(GeometryPointString outerSlipPlane)
        {
            outerSlipPlane.Points.Clear();

            GeometryPoint dikeTopRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeToePolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint bottomDitchPolderSide = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchPolderSide);
            GeometryPoint bottomDitchDikeSide = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchDikeSide);
            GeometryPoint surfaceLevelInside = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside);

            var exitPoint = new GeometryPoint();
            bool isDitchPresent = IsDitchPresent();

            // Entry point (1out)
            GeometryPoint entryPoint = dikeTopRiver;

            //  Height h is calculated at DikeToeAtPolder
            var h = dikeToePolder.Z - GetTopPleistoceen(dikeToePolder.X);

            // If h is nul, the height between surface line and bottom boundary is used
            if (h < Constants.CAlmostZero)
            {
                h = dikeToePolder.Z - profile2D.Geometry.Bottom;
            }
            // Exit point (9out) and level of the bottom part (6out to 8out)
            double bottom;
            if (isDitchPresent)
            {
                if (dikeToePolder.Z - 2 * h / 3 > bottomDitchDikeSide.Z || dikeToePolder.Z - 2 * h / 3 > bottomDitchPolderSide.Z)
                {
                    exitPoint = bottomDitchDikeSide; // point O
                    bottom = bottomDitchDikeSide.Z; // point O
                }
                else
                {
                    exitPoint = bottomDitchPolderSide; // point P
                    bottom = dikeToePolder.Z - 2 * h / 3;
                }
            }
            else
            {
                exitPoint.X = Math.Min(dikeToePolder.X + 20.0, surfaceLevelInside.X - 0.01);
                exitPoint.Z = surfaceLine.Geometry.GetZAtX(exitPoint.X);
                bottom = Math.Min(exitPoint.Z, dikeToePolder.Z - 2 * h / 3);
            }

            // determine corner point: intersection point between line 1 (vertical) and the horizontal line
            var cornerPoint = new GeometryPoint();

            var vertLength = entryPoint.Z - bottom;
            var vertSegmentLength = vertLength / 4;
            cornerPoint.X = entryPoint.X;
            cornerPoint.Z = entryPoint.Z - 4 * vertSegmentLength;

            outerSlipPlane.Points.Add(entryPoint);

            // Determine the intersection of the final slope and the vertical part
            var intersection = new GeometryPoint();
            intersection.Z = exitPoint.Z - (exitPoint.X - entryPoint.X);
            bool isDeepPleistocene = intersection.Z >= bottom;

            // Normal configuration: Outer slip plane has an horizontal bottom part
            if (!isDeepPleistocene)
            {
                // The 4 points on the first (vertical) line    
                for (int i = 1; i < 4; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = entryPoint.X;
                    nextPoint.Z = entryPoint.Z - i * vertSegmentLength;
                    outerSlipPlane.Points.Add(nextPoint);
                }

                // Point 8out
                var point8out = new GeometryPoint();
                if (Math.Abs(exitPoint.Z - bottom) < Constants.CGeoAccu)
                {
                    point8out.X = exitPoint.X - (exitPoint.X - entryPoint.X) / 4;
                }
                else
                {
                    point8out.X = exitPoint.X - (exitPoint.Z - bottom); // last part has a slope 1:1
                }
                point8out.Z = bottom;

                // Point 6out
                var point6out = new GeometryPoint();

                var horLength = point8out.X - entryPoint.X;
                var horSegmLength = horLength / 3;
                point6out.X = entryPoint.X + horSegmLength;
                point6out.Z = bottom;

                // Point 5out
                var point5out = new GeometryPoint();
                point5out.X = entryPoint.X + horSegmLength / 2;
                point5out.Z = bottom + vertSegmentLength / 2;

                // Point 7out
                var point7out = new GeometryPoint();
                point7out.X = point6out.X + horSegmLength;
                point7out.Z = bottom;

                outerSlipPlane.Points.Add(point5out);
                outerSlipPlane.Points.Add(point6out);
                outerSlipPlane.Points.Add(point7out);
                outerSlipPlane.Points.Add(point8out);
                outerSlipPlane.Points.Add(exitPoint);
            }
            // Other configuration (i.e. deep Pleistocene): Outer slip plane has no horizontal bottom part, but only a first vertical part and an oblique part with slope 1:1
            else
            {
                // The 5 points on the first vertical line (= points 1out to 5out included)
                for (int i = 1; i < 5; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = entryPoint.X;
                    nextPoint.Z = entryPoint.Z - i * (entryPoint.Z - intersection.Z) / 4;
                    outerSlipPlane.Points.Add(nextPoint);
                }

                double vertSegmentLengthForLastSlope = (exitPoint.Z - intersection.Z) / 4;
                // The 4 last points (i.e. 6out to 9out)
                for (int i = 1; i < 5; i++)
                {
                    var nextPoint = new GeometryPoint();
                    nextPoint.X = entryPoint.X + i * vertSegmentLengthForLastSlope;
                    nextPoint.Z = intersection.Z + i * vertSegmentLengthForLastSlope;
                    outerSlipPlane.Points.Add(nextPoint);
                }
            }

            outerSlipPlane = EnsureThatPointsOfTheSlipPlaneAreNotAboveSurfaceLine(outerSlipPlane);

            DataEventPublisher.DataListModified(outerSlipPlane.Points);
        }

        /// <summary>
        /// If points of the Spencer slip plane are above the surface line, they are shifted 10 cm below the surface line
        /// </summary>
        /// <param name="slipPlane">Slip plane</param>
        /// <returns>Updated slip plane</returns>
        private GeometryPointString EnsureThatPointsOfTheSlipPlaneAreNotAboveSurfaceLine(GeometryPointString slipPlane)
        {
            foreach (var points in slipPlane.Points)
            {
                double surfaceLevel = surfaceLine.Geometry.GetZAtX(points.X);
                if (points.Z - surfaceLevel > Constants.CGeoAccu)
                {
                    points.Z = surfaceLevel - 0.1;
                }
            }
            return slipPlane;
        }

        private void UpdateSpencerGA(StabilityModel model)
        {
            GeometryPoint dikeTopRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeTopPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            double xMiddleDike = (dikeTopPolder.X + dikeTopRiver.X) / 2;

            if (!double.IsNaN(GetTopPleistoceen(xMiddleDike)))
            {
                if (model.IsHighSlipPlanes)
                {
                    CreateInwardSpencerPlanesHigh(model);
                }
                else
                {
                    CreateInwardSpencerPlanesLow(model);
                }
            }
        }

        private void CreateInwardSpencerPlanesHigh(StabilityModel model)
        {
            CreateSlipPlanes(model);

            CreateInwardSpencerInnerSlipPlaneHigh(model.SlipPlanes[0]);
            CreateInwardSpencerOuterSlipPlaneHigh(model.SlipPlanes[1]);
        }

        private void CreateSlipPlanes(StabilityModel model)
        {
            while (model.SlipPlanes.Count < 2)
            {
                model.CreateSlipPlane(new List<GeometryPoint>(), true);
            }
        }

        private void CreateInwardSpencerPlanesLow(StabilityModel model)
        {
            CreateSlipPlanes(model);

            CreateInwardSpencerInnerSlipPlaneLow(model.SlipPlanes[0]);
            CreateInwardSpencerOuterSlipPlaneLow(model.SlipPlanes[1]);
        }

        /// <summary>
        /// Finds the intersection point of a given horizontal line and the river side talud.
        /// </summary>
        /// <param name="surfaceLine">The evaluated surfaceline.</param>
        /// <param name="obliqueLine">The oblique line.</param>
        /// <returns>The intersection point, or null in case no intersection was found.</returns>
        private static GeometryPoint DetermineIntersectionBetweenDikeSlopeAtRiverAndAnObliqueLine(SurfaceLine2 surfaceLine, Line obliqueLine)
        {
            var list = surfaceLine.Geometry.Points.Where(point =>
                                                         point.X >= surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside).X &&
                                                         point.X <= surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X).ToList();

            list.Insert(0, surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside));
            list.Add(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver));
            list.Sort();

            for (int i = 1; i < list.Count; i++)
            {
                var surfaceLineSegment = new Line();
                surfaceLineSegment.SetBeginAndEndPoints(
                    new GeometryPoint(list[i - 1].X, 0, list[i - 1].Z),
                    new GeometryPoint(list[i].X, 0, list[i].Z));

                var intersectPoint = surfaceLineSegment.GetIntersectPointXZ(obliqueLine);
                if (intersectPoint != null)
                {
                    return GeometryPoint.CreateNewXZPoint(intersectPoint.X, intersectPoint.Z);
                }
            }

            return null;
        }

        private void ResetSlipPlanes(StabilityModel stabilityModel)
        {
            stabilityModel.SlipPlanes.Clear();
            stabilityModel.SlipPlanes.Add(new SlipPlane
                {
                    Name = LocalizationManager.GetTranslatedText(typeof(StabilityDikeLocationInfo), "SlipPlane")
                }); // One empty slip plane is expected to be present by default
        }

        /// <summary>
        /// Search the X coordinate at which the uplift ends, starting the search at Dike Toe At Polder with an increment of 50 cm.
        /// If uplift does not occurs, this procedure returns the X coordinate at Surface Level Inside.
        /// </summary>
        /// <returns>X coordinate at which uplift ends</returns>
        private double SearchXcoordinateOfEndOfUpliftAtPolderSide()
        {
            const double searchUpliftIncrement = 0.5;
            double xStart = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
            double xLast = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X;
            bool isUplift = false;

            var upliftCalculator = new UpliftCalculator();
            upliftCalculator.SoilProfile2D = profile2D;
            upliftCalculator.Waternet = waternet;

            // First search at which X coordinate uplift starts
            while ((!isUplift) && ((xStart + searchUpliftIncrement) < xLast))
            {
                xStart += searchUpliftIncrement;
                double upliftPotential = upliftCalculator.DetermineUpliftPotentialAtXAtTopOfHighestAquiferFromA2DProfile(xStart);
                isUplift = (upliftPotential < 1.0);
            }
            // If no uplift occurs, return the X coordinate at Surface Level Inside.
            if (isUplift == false)
            {
                return xLast;
            }
            // If uplift occurs, search the X coordinate at which uplift ends
            while ((isUplift) && ((xStart + searchUpliftIncrement) < xLast))
            {
                xStart += searchUpliftIncrement;
                double upliftPotential = upliftCalculator.DetermineUpliftPotentialAtXAtTopOfHighestAquiferFromA2DProfile(xStart);
                isUplift = (upliftPotential < 1.0);
            }
            return xStart;
        }
    }
}