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

namespace Deltares.Geotechnics.WaternetCreator
{
    public class ExpertKnowledgePLLinesCreator : PlLinesCreator
    {
        public void CreatePL1() { }
    }

    public class RingtoetsWTI2017PLLinesCreator
    {
        public void CreatePL1() { }
    }

    public class DupuitPLLinesCreator
    {
        public void CreatePL1() { }
    }

    public class PlLinesCreator
    {
        private const double cUpliftFactorEquilibrium = 1.0;
        //TODO : ROB constants?
        private const double cEpsilon = 1.0e-3;
        private const double cOffsetPhreaticLineBelowSurface = 0.01;
        private const int cRoundDecimals = 6;
        private double HeadAtDikeToePolderPlus20M;
        private double HeadAtEndOfGeometry;
        private GeometryPointString Pl2;
        private GeometryPointString phreaticLine;
        //private double waterLevelRiverHigh = 0;

        /*public PhreaticAdaptionType? NWOPhreaticAdaption { get; set; }*/

        // Output
        private double pl3HeadAdjusted;
        private double pl3LocationXMinUplift;
        private double pl3MinUplift;
        private double pl4HeadAdjusted;
        private double pl4LocationXMinUplift;
        private double pl4MinUplift;
        private Dictionary<double, double> upliftValuesPl3;
        private Dictionary<double, double> upliftValuesPl4;
        private double waterLevelPolder = 0;

        /// <summary>
        ///     Constructor
        /// </summary>
        public PlLinesCreator()
        {
            InitPLLinesCreator();
        }

        public bool IsUseOvenDryUnitWeight { get; set; }
        public SoilGeometryType SoilGeometryType { get; set; }
        public string SoilGeometry2DName { get; set; }
        public Soil DikeEmbankmentMaterial { get; set; }

        //public SoilProfile1D SoilProfile { get; set; }
        //public SoilProfile2D SoilProfile2D { get; set; }
        public double? HeadInPLLine2 { get; set; }
        public double GaugeMissVal { get; set; }

        public double Pl3HeadAdjusted
        {
            get
            {
                return pl3HeadAdjusted;
            }
        }

        public double Pl3LocationXMinUplift
        {
            get
            {
                return pl3LocationXMinUplift;
            }
        }

        public double Pl4HeadAdjusted
        {
            get
            {
                return pl4HeadAdjusted;
            }
        }

        public double Pl4LocationXMinUplift
        {
            get
            {
                return pl4LocationXMinUplift;
            }
        }

        public bool Inwards { get; set; }

        /// <summary>
        /// Create all PL-lines
        /// </summary>
        /// <param name="location">Class containing all hydraulic data's needed for PL-line creation</param>
        /// <returns>PL1 (phreatic line), PL2, PL3 and PL4</returns>
        public Dictionary<PLLineType, GeometryPointString> CreateAllPLLines(Location location)
        {
            var plLines = new Dictionary<PLLineType, GeometryPointString>();
            switch (location.PlLineCreationMethod)
            {
                case PlLineCreationMethod.ExpertKnowledgeLinearInDike:
                case PlLineCreationMethod.ExpertKnowledgeRrd:
                case PlLineCreationMethod.RingtoetsWti2017:
                    plLines = CreateAllPLLinesWithExpertKnowledge(location);
                    break;
            }

            if (!LineHelper.IsXAscending(plLines[PLLineType.PL1].Points))
            {
                string message = "PLLine 1 not an X-ascending polyline";
                throw new System.Exception(message);
            }

            return plLines;
        }

        /// <summary>Create 4 PL-lines for the waternet creator of Ringtoets WTI 2017.
        ///  </summary>
        /// <param name="plLineType">Enter the PL type you want to create</param>
        /// <param name="location">Class containing all hydraulic data needed for the creation of PL-lines</param>
        /// <returns>All the points of the desired PL-line></returns>
        public GeometryPointString CreatePLLineByExpertKnowledge(PLLineType plLineType, Location location)
        {
            GeometryPointString plLine = null;

            switch (plLineType)
            {
                case PLLineType.PL1:
                    plLine = CreatePLLine1ByExpertKnowledge(location);
                    break;
                case PLLineType.PL2:
                    plLine = CreatePLLine2ByExpertKnowledge(location);
                    break;
                case PLLineType.PL3:
                    plLine = CreatePLLine3ByExpertKnowledge(location);
                    break;
                case PLLineType.PL4:
                    plLine = CreatePLLine4ByExpertKnowledge(location);
                    break;
            }

            if (plLine != null)
            {
                LineHelper.DeleteCoinsidingPoints(plLine.Points, 0.0001);
            }
            return plLine;
        }

        /// <summary>
        ///     Create PL2 (is pl line for penetration)
        /// </summary>
        /// <returns></returns>
        public GeometryPointString CreatePLLine2ByExpertKnowledge(Location location)
        {
            GeometryPointString plLine = null;
            if ((location.PenetrationLength.AlmostEquals(0.0)) || Double.IsNaN(location.HeadInPLLine2Outwards))
            {
                // No penetration, or no Head Pl2 defined, so empty pl-line will be returned
                plLine = new GeometryPointString();
            }
            else
            {
                plLine = CreatePLLine2ByExpertKnowledgeFor1DGeometry(location);
            }

            return plLine;
        }

        /// <summary>
        /// Create PL-line 3 (= head in the Pleistocene)
        /// </summary>
        /// <returns>PL-line 3</returns>
        public GeometryPointString CreatePLLine3ByExpertKnowledge(Location location)
        {
            //TODO ROB: check of HeadInPLLine3 != double.nan ?
            if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017)
            {
                switch (location.DikeSoilScenario)
                {
                    case DikeSoilScenario.ClayDikeOnClay:
                    case DikeSoilScenario.SandDikeOnClay:
                        return CreatePLLine3Or4ForRingtoetsWti2017ForDikesOnClay(location.HeadInPLLine3, PLLineType.PL3,
                                                                                 location);
                    case DikeSoilScenario.ClayDikeOnSand:
                        return createPLLine3ClayDikeOnSandForRingtoetsWti2017Method(location);
                }
            }
            else
            {
                switch (location.DikeSoilScenario)
                {
                    case DikeSoilScenario.ClayDikeOnClay:
                        return CreatePLLine3Or4ByExpertKnowledge(location.HeadInPLLine3, PLLineType.PL3, location);
                    case DikeSoilScenario.SandDikeOnClay:
                        return CreatePLLine3Or4ByExpertKnowledge(location.HeadInPLLine3, PLLineType.PL3, location);
                    case DikeSoilScenario.ClayDikeOnSand:
                        return createPLLine3ClayDikeOnSandForRingtoetsWti2017Method(location);
                }
            }
            return null;
        }

        /// <summary>
        /// Create PL-line 4 (= head in the in-between aquifer)
        /// </summary>
        /// <returns>PL-line 4</returns>
        public GeometryPointString CreatePLLine4ByExpertKnowledge(Location location)
        {
            if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017)
            {
                switch (location.DikeSoilScenario)
                {
                    case DikeSoilScenario.ClayDikeOnClay:
                    case DikeSoilScenario.SandDikeOnClay:
                        return CreatePLLine3Or4ForRingtoetsWti2017ForDikesOnClay(location.HeadInPLLine4, PLLineType.PL4, location);
                }
            }
            else
            {
                return CreatePLLine3Or4ByExpertKnowledge(location.HeadInPLLine3, PLLineType.PL4, location);
            }
            return null;
        }

        /// <summary>
        /// Create PL-line 3 (= head in Pleistocene) for "Clay dike on sand" (used only for Ringtoets WTI 2017 method)
        /// </summary>
        /// <param name="location"></param>
        /// <returns></returns>
        public GeometryPointString CreatePLLine3ClayDikeOnSand(Location location)
        {
            double headValue = location.WaterLevelRiver;
            var plDiketopRiver =
                new GeometryPoint(location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X, 0, headValue);
            var plDiketoepolder =
                new GeometryPoint(location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, 0,
                                  location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).Z);
            var plLine = new GeometryPointString();
            double firstX = location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X;
            plLine.Points.Add(new GeometryPoint(firstX, 0, headValue));
            plLine.Points.Add(plDiketopRiver);
            plLine.Points.Add(plDiketoepolder);
            double gradient = (plDiketopRiver.X - plDiketoepolder.X) / (plDiketopRiver.Z - plDiketoepolder.Z);
            double xDistanceDiketoeToPolderpeilIntersection = (location.WaterLevelPolder - plDiketoepolder.Z) * gradient;
            plLine.Points.Add(new GeometryPoint(plDiketoepolder.X + xDistanceDiketoeToPolderpeilIntersection, 0,
                                                location.WaterLevelPolder));

            double lastPointX = location.Surfaceline.Geometry.Points.Max(p => p.X);
            plLine.Points.Add(new GeometryPoint(lastPointX, 0, location.WaterLevelPolder));

            plLine.CondensePoints();
            return plLine;
        }

        /// <summary>
        ///     get uplift factors for x, regarding pl4
        /// </summary>
        /// <returns></returns>
        public Dictionary<double, double> GetUpliftOutputValuesPl4()
        {
            return upliftValuesPl4;
        }

        /// <summary>
        ///     get uplift factors for x, regarding pl3
        /// </summary>
        /// <returns></returns>
        public Dictionary<double, double> GetUpliftOutputValuesPl3()
        {
            return upliftValuesPl3;
        }

        public bool IsNonWaterRetainingObjectPoint(GeometryPoint point, SurfaceLine2 surfaceline)
        {
            return
                ((point.LocationEquals(
                    surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1))) ||
                 (point.LocationEquals(
                     surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint2))) ||
                 (point.LocationEquals(
                     surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint3))) ||
                 (point.LocationEquals(
                     surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4)))
                );
        }

        /// <summary>
        ///     Intersection between two line segments:
        ///     - Horizontal waterlevel
        ///     - Surface line (until DikeTopAtRiver)
        /// </summary>
        /// <param name="waterLevelRiver"></param>
        /// <param name="surfaceline"> </param>
        /// <returns></returns>
        public GeometryPoint DetermineIntersectionBetweenWaterLevelAndDike(double waterLevelRiver,
                                                                           SurfaceLine2 surfaceline)
        {
            GeometryPoint pointAtLevel = DetermineIntersectionBetweenTaludRiverSideAndWaterLevel(waterLevelRiver,
                                                                                                 surfaceline);
            //.DetermineIntersectionWithLevel(waterLevelRiver);
            if (pointAtLevel != null)
            {
                return new GeometryPoint(pointAtLevel.X, 0, pointAtLevel.Z);
            }
            return null;
        }

        /// <summary>
        ///     Determines the intersection between talud river side and water level.
        ///     Search starts from DikeTopAtRiver to SurfaceLevelOutside
        /// </summary>
        /// <param name="level">The water level or z value of the horizontal line</param>
        /// <param name="surfaceline"> </param>
        /// <returns>An intersection point if any null otherswise</returns>
        public GeometryPoint DetermineIntersectionBetweenTaludRiverSideAndWaterLevel(double level, SurfaceLine2 surfaceline)
        {
            var waterlevelLine = new Line();
            double xStart = surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside).X;
            double xEnd = surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X;
            waterlevelLine.CreateHorizontalZLine(xStart, xEnd, level);

            ThrowWhenWaterLevelAboveDike(level, surfaceline);

            var relevantPart = new GeometryPointString();
            relevantPart.Points.AddRange(surfaceline.Geometry.Points.Where(p => !double.IsNaN(p.X)));

            double x = relevantPart.GetXAtZ(level);

            if (x <= surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X + cEpsilon)
            {
                return new GeometryPoint(x, 0, level);
            }
            else
            {
                return null;
            }
        }

        public GeometryPoint DetermineIntersectionBetweenPl2Or1AndFinalSlope100OfPl3(GeometryPoint firstPointOfObliqueLine, double slopeOfObliqueLine, Location location)
        {
            // Creation of the oblique line with a slope of slopeOfObliqueLine (1:100)
            var obliqueLineWithSlope100 = new Line();
            const double deltaX = 10000;
            double signSlope;
            if (CalculateHeadPl2AtX(location, firstPointOfObliqueLine.X) < firstPointOfObliqueLine.Z)
            {
                signSlope = -1;
            }
            else
            {
                signSlope = 1;
            }
            double xFirstPointOfObliqueLine = firstPointOfObliqueLine.X - deltaX;
            double zFirstPointOfObliqueLine = firstPointOfObliqueLine.Z + signSlope * (xFirstPointOfObliqueLine - firstPointOfObliqueLine.X) / slopeOfObliqueLine;
            double xLastPointOfObliqueLine = firstPointOfObliqueLine.X + deltaX;
            double zLastPointOfObliqueLine = firstPointOfObliqueLine.Z + signSlope * (xLastPointOfObliqueLine - firstPointOfObliqueLine.X) / slopeOfObliqueLine;

            obliqueLineWithSlope100.SetBeginAndEndPoints(new GeometryPoint(xFirstPointOfObliqueLine, 0, zFirstPointOfObliqueLine),
                new GeometryPoint(xLastPointOfObliqueLine, 0, zLastPointOfObliqueLine));

            // Creation of the relevant part of PL2 (or PL1 if PL2 is not present)
            var segmentPl2OrPl1 = new Line();
            if (Pl2 != null)
            {
                double x1 = firstPointOfObliqueLine.X - deltaX;
                double x2 = firstPointOfObliqueLine.X + deltaX;
                double z1 = CalculateHeadPl2AtX(location, x1);
                double z2 = CalculateHeadPl2AtX(location, x2);
                segmentPl2OrPl1.SetBeginAndEndPoints(new GeometryPoint(x1, 0, z1), new GeometryPoint(x2, 0, z2));
                GeometryPoint intersectPoint = LineHelper.GetIntersectPoint(segmentPl2OrPl1, obliqueLineWithSlope100, PointType.XZ);
                if (!double.IsNaN(intersectPoint.X))
                {
                    return new GeometryPoint(intersectPoint.X, 0, intersectPoint.Z);
                }
            }
            else
            {
                GeometryPointString relevantPartOfPl1 = phreaticLine.GetPart(firstPointOfObliqueLine.X, firstPointOfObliqueLine.X + deltaX);
                for (int index = 0; index < relevantPartOfPl1.Points.Count - 1; index++)
                {
                    segmentPl2OrPl1 = new Line(relevantPartOfPl1[index], relevantPartOfPl1[index + 1]);
                    GeometryPoint intersectPoint = LineHelper.GetIntersectPoint(segmentPl2OrPl1, obliqueLineWithSlope100, PointType.XZ);
                    if (!double.IsNaN(intersectPoint.X))
                    {
                        return new GeometryPoint(intersectPoint.X, 0, intersectPoint.Z);
                    }
                }
            }
            return null;
        }

        public static bool Equal(double aValue1, double aValue2)
        {
            if (Math.Abs(aValue1 - aValue2) < cEpsilon)
            {
                return true;
            }

            return false;
        }

        /// <summary>
        /// Create PL1 (= phreatic level)
        /// For methods Expert Knowledge, refer to the user manual of DAM. 
        /// For method Ringtoets WTI 2017 , refer to the memo of Alexander van Duinen in:
        /// https://repos.deltares.nl/repos/delftgeosystems/delftgeosystems/trunk/documents/Geotechnics
        /// </summary>
        /// <returns>PL1 (= phreatic level)</returns>
        public GeometryPointString CreatePLLine1ByExpertKnowledge(Location location)
        {
            /*this validation is carried out within location*/
            //ValidateSurfaceLine(location.Surfaceline);
            //ValidateRequiredCharacteristicPoints(location.Surfaceline);
            //ThrowWhenWaterLevelPolderAboveDikeTopAtPolder(location);
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeTopAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);

            //Create Phreatic line and add polderwater level            
            double startXCoordinate =
                location.Surfaceline.Geometry.Points.Where(p => !double.IsNaN(p.X)).OrderBy(p => p.X).First().X;

            var line = new GeometryPointString();
            double waterLevel;
            if (Inwards)
            {
                waterLevel = location.WaterLevelRiver;
            }
            else
            {
                waterLevel = location.WaterLevelRiverLow;
            }
            // double waterLevelRiver = 0;
            /*if (double.IsNaN(location.WaterLevelRiver))
            {
               location.WaterLevelRiver = !double.IsNaN(location.WaterLevelRiverLow) ? location.WaterLevelRiverLow : location.WaterLevelRiverAverage;
            }*/

            // Water level at river side is supposed to be at level of SurfaceLevelOutside when waterlevel is below SurfaceLevelOutside (is the bottom of the river)
            if (!double.IsNaN(EnsureWaterLevelIsAboveRiverBottom(waterLevel, location.Surfaceline)))
            {
                waterLevel = EnsureWaterLevelIsAboveRiverBottom(waterLevel, location.Surfaceline);
                //location.WaterLevelRiver = waterLevelRiver;
            }

            line.Points.Add(new GeometryPoint(startXCoordinate, 0, waterLevel));

            GeometryPoint intersectionPhreaticDike = null;

            // Determine intersection between the Phreatic Level and the berm (or dike) at river side 
            intersectionPhreaticDike = DetermineIntersectionBetweenWaterLevelAndDike(waterLevel, location.Surfaceline);
            if (intersectionPhreaticDike == null)
            {
                string message = "DetermineIntersectionBetweenWaterLevelAndDike failed";
                throw new System.Exception(message);
            }

            // Add intersection point (i.e. point A of figure 19 handleiding DAM 1.0 - part C) to phreatic line list if not already in list.
            if (!line.Points.Contains(intersectionPhreaticDike))
            {
                line.Points.Add(intersectionPhreaticDike);
            }

            // Complete phreatic line with points below dike crest (i.e. dike top at river and dike top at polder)                                                          
            if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017)
            {
                CreatePhreaticLineSegmentsInsideDikeForRingtoetsWti2017(line, location);
            }
            else
            {
                if (Inwards)
                {
                    CreatePhreaticLineSegmentsInsideDikeForHighRiverLevel(line, location);
                }
                else
                {
                    // intersectionPhreaticDike.Z is either low riverlevel or ToeAtDike.Z, whichever is higher
                    // Or in human language: low riverlevel should always be a least at surfacelevel
                    CreatePhreaticLineSegmentsInsideDikeForLowRiverLevel(line, intersectionPhreaticDike.Z, location);
                }
            }

            if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017)
            {
                CreatePhreaticLineSegmentsInShoulderAndPolderForRingtoetsWti2017(line, location);
            }
            else
            {
                CreatePhreaticLineSegmentsInShoulderAndPolder(line, location);
            }

            // Check if phreatic line is above 
            line.CondensePoints();
            ValidatePhreaticAboveWaterLevelPolder(line, location);
            // Check if phreatic line is below the surface line
            ValidatePhreaticLineBelowDike(line, location);

            // For Ringtoets WTI 2017 method, both points at dike top have a minimum value
            if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017)
            {
                line.CondensePoints();
                line.SortPointsByXAscending();
                // First dike top at river is created for all scenario's
                double interpolatedLevelDikeTopRiver = line.GetZAtX(dikeTopAtRiver.X);
                double waterLevelBelowDikeTopAtRiver = Math.Max(interpolatedLevelDikeTopRiver, location.MinimumLevelPhreaticLineAtDikeTopRiver);
                line.Points.Add(new GeometryPoint(dikeTopAtRiver.X, 0, waterLevelBelowDikeTopAtRiver));
                line.SortPointsByXAscending();

                // Then dike top at polder is created only for sand dikes (for clay dikes it was already created
                // in procedure CreatePhreaticLineSegmentsInsideDikeForRingtoetsWti2017)
                if ((location.DikeSoilScenario == DikeSoilScenario.SandDikeOnClay || location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand)
                    && location.UseDefaultOffsets)
                {
                    double interpolatedLevelDikeTopPolder = line.GetZAtX(dikeTopAtPolder.X);
                    double waterLevelBelowDikeTopAtPolder = Math.Max(interpolatedLevelDikeTopPolder, location.MinimumLevelPhreaticLineAtDikeTopPolder);
                    line.Points.Add(new GeometryPoint(dikeTopAtPolder.X, 0, waterLevelBelowDikeTopAtPolder));
                    line.SortPointsByXAscending();
                }
                ValidatePhreaticLineBelowDike(line, location);
            }

            // currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
            AdaptPL1ForNonWaterRetainingObject(line, location);

            // remove all double points
            line.CondensePoints();

            phreaticLine = line;
            return line;
        }

        /// <summary>
        /// Make sure that river level is above the bottom of the river
        /// </summary>
        /// <param name="riverLevel"></param>
        /// <param name="surfaceline"> </param>
        /// <returns>bottom of the river if it is above riverLevel, else riverLevel</returns>
        public double EnsureWaterLevelIsAboveRiverBottom(double riverLevel, SurfaceLine2 surfaceline)
        {
            // Waterlevel is supposed to be at level of SurfaceLevelOutside when waterlevel
            // is below SurfaceLevelOutside (is the bottom of the river)
            if (!double.IsNaN(riverLevel))
            {
                double z = surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside).Z;
                riverLevel = Math.Max(z, riverLevel);
            }
            return riverLevel;
        }

        /// <summary>
        ///     If shoulder is present then the toe of the shoulder is returned
        ///     Else the toe of the dike (ShoulderSaseInside) is returned
        /// </summary>
        /// <param name="surfaceline"> </param>
        /// <returns>toe of the dike</returns>
        public GeometryPoint GetDikeToeInward(SurfaceLine2 surfaceline)
        {
            if (surfaceline.IsDefined(CharacteristicPointType.ShoulderBaseInside))
            {
                return surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside);
            }
            return surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
        }

        public double DetermineZValueAtX(double x, double x1, double x2, double z1, double z2)
        {
            double fraction = (z2 - z1) / (x2 - x1);
            return fraction * (x - x1) + z1;
        }

        /// <summary>
        ///     Initialize
        /// </summary>
        private void InitPLLinesCreator()
        {
            IsUseOvenDryUnitWeight = false;
        }

        /// <summary>
        /// Create the PL-lines for ExpertKnowledge method (used in DAM) and Ringtoest WTI 2017 method
        /// </summary>
        /// <param name="location">Class containing all hydraulic data's needed for PL-line creation</param>
        /// <returns>PL1 (phreatic line), PL2, PL3 and PL4</returns>
        private Dictionary<PLLineType, GeometryPointString> CreateAllPLLinesWithExpertKnowledge(Location location)
        {
            var plLines = new Dictionary<PLLineType, GeometryPointString>();
            foreach (PLLineType plLineType in Enum.GetValues(typeof(PLLineType)))
            {
                bool isPL1LineDefinedForLocation = (location != null) && (location.LocalXzpl1Line != null) &&
                                                   (location.LocalXzpl1Line.Points.Count > 1);
                if ((plLineType == PLLineType.PL1) && isPL1LineDefinedForLocation)
                {
                    GeometryPointString plLine = plLines[plLineType];
                }
                else
                {
                    plLines[plLineType] = CreatePLLineByExpertKnowledge(plLineType, location);
                }
            }

            return plLines;
        }

        /// <summary>
        /// Get soil profile 1D below <param name="xCoordinate"></param> using <param name="soilProfile2D"></param>.
        /// </summary>
        /// <param name="soilProfile2D">Soil profile 2D</param>
        /// <param name="xCoordinate">X coordinate</param>
        /// <returns>The created soil profile 1D </returns>
        private SoilProfile1D GetSoilProfileBelowPoint(SoilProfile2D soilProfile2D, double xCoordinate)
        {
            return soilProfile2D.GetSoilProfile1D(xCoordinate);
        }

        /// <summary>
        /// Create PL2 (is PL-line for penetration) for 1D-geometry: the soil geometry at Dike Top at Polder is used and PL2 is created
        /// only if an aquitard above an aquifer is found. If created, PL-line 2 is a straight line from left to right.
        /// </summary>
        /// <param name="location">Class containing all the hydraulic data's</param>
        /// <returns>PL-line 2</returns>
        private GeometryPointString CreatePLLine2ByExpertKnowledgeFor1DGeometry(Location location)
        {
            SoilProfile1D actualSoilProfile = location.SoilProfile2D != null
                                                  ? GetSoilProfileBelowPoint(location.SoilProfile2D,
                                                                             location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).X)
                                                  : location.SoilProfile1D;

            var plLine = new GeometryPointString();

            if (actualSoilProfile != null)
            {
                IList<SoilLayer1D> aquiferLayers = actualSoilProfile.GetAquiferLayers();

                // ROB: Add validation method to check this
                if (aquiferLayers.Count == 0)
                {
                    string message = "Soil profile contains no aquifer layers.";
                    throw new System.Exception(message);
                }

                if (location.PenetrationLength > 0 && aquiferLayers.Count > 0)
                {
                    IList<SoilLayer1D> infiltrationLayers = (from SoilLayer1D layer in actualSoilProfile.Layers
                                                             where
                                                                 (actualSoilProfile.InBetweenAquiferLayer == null ||
                                                                  layer.TopLevel <
                                                                  actualSoilProfile.InBetweenAquiferLayer.TopLevel) &&
                                                                 layer.TopLevel >
                                                                 actualSoilProfile.BottomAquiferLayer.TopLevel
                                                             select layer).ToList();

                    if (infiltrationLayers.Count > 0)
                    {
                        double separationLevel = actualSoilProfile.BottomAquiferLayer.TopLevel +
                                                 location.PenetrationLength;

                        if (separationLevel <= infiltrationLayers.First().TopLevel)
                        {
                            plLine.Points.Add(
                                new GeometryPoint(location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X, 0,
                                                  location.CorrectedHeadInPLLine2Outwards));
                            plLine.Points.Add(
                                new GeometryPoint(location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X, 0,
                                                  location.CorrectedHeadInPLLine2Inwards));
                        }
                    }
                }
            }
            Pl2 = plLine;
            return plLine;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="location"></param>
        /// <returns></returns>
        private GeometryPointString createPLLine3ClayDikeOnSandForRingtoetsWti2017Method(Location location)
        {
            var plLine = new GeometryPointString();
            if (!location.Surfaceline.HasAnnotation(CharacteristicPointType.DikeToeAtRiver))
            {
                throw new System.Exception("Characteristic point \"dike toe at river\" is not defined.");
            }
            var outerBoundary = new GeometryPoint();
            outerBoundary.X = location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X;
            outerBoundary.Z = location.WaterLevelRiver;
            plLine.Points.Add(outerBoundary);

            var outertoplevel = new GeometryPoint();
            outertoplevel.X = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X;
            outertoplevel.Z = location.WaterLevelRiver;
            plLine.Points.Add(outertoplevel);

            var innerToe = new GeometryPoint();
            innerToe.X = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
            innerToe.Z = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).Z;
            plLine.Points.Add(innerToe);

            var innerBoundary = new GeometryPoint();
            innerBoundary.X = location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X;
            innerBoundary.Z = location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).Z;
            plLine.Points.Add(innerBoundary);

            plLine.CondensePoints();
            return plLine;
        }

        private GeometryPointString CreatePLLine3Or4ByExpertKnowledge(double headValue, PLLineType plLineType, Location location)
        {
            var plLine = new GeometryPointString();
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeToeAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
            GeometryPoint dikeTopAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint dikeToeAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);

            // Soilprofile is used to check if there is really an aquifer below toe of dike at river.
            // The assumption is made that if there is an aquifer at the riverside and at the the polderside,
            // that there is a connection between these aquifers.
            // In the uplift calculation there will also be a check on the existence of an aquifer.
            // 
            if (!location.Surfaceline.HasAnnotation(CharacteristicPointType.DikeToeAtRiver))
            {
                throw new System.Exception("Characteristic point \"dike toe at river\" is not defined.");
            }
            SoilProfile1D actualSoilProfile = location.SoilProfile2D != null
                                                  ? GetSoilProfileBelowPoint(location.SoilProfile2D, dikeToeAtRiver.X) : location.SoilProfile1D;
            if (actualSoilProfile == null)
            {
                throw new System.Exception("No soil profile is present.");
            }

            SoilLayer1D relevantAquiferLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile, location);

            if (relevantAquiferLayer != null && !double.IsNaN(location.LeakageLengthInwardsPl3))
            {
                // initialize
                double pl2level = 0;
                double pl2gradient = 0;
                double firstX = location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X;
                double lastX = location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X;
                double headInPlLine2Outwards = location.CorrectedHeadInPLLine2Outwards;
                if (double.IsNaN(headInPlLine2Outwards))
                {
                    headInPlLine2Outwards = location.WaterLevelPolder;
                    pl2level = location.WaterLevelPolder;
                }
                else
                {
                    //gradient is calculated from dike toe out to dike toe inwards
                    pl2gradient = (headInPlLine2Outwards - location.CorrectedHeadInPLLine2Inwards) / (firstX - lastX);
                    pl2level = headInPlLine2Outwards + pl2gradient * (dikeTopAtRiver.X - firstX);
                }
                double waterLevelRiverAverage = !double.IsNaN(location.WaterLevelRiverAverage)
                                                    ? location.WaterLevelRiverAverage
                                                    : headInPlLine2Outwards;
                double headPlDikeTopRiver = 0;
                double LeakageLengthOutward = 0;
                double LeakageLengthInward = 0;
                if (plLineType == PLLineType.PL3)
                {
                    LeakageLengthOutward = location.LeakageLengthOutwardsPl3;
                    LeakageLengthInward = location.LeakageLengthInwardsPl3;
                }
                else
                {
                    LeakageLengthOutward = location.LeakageLengthOutwardsPl4;
                    LeakageLengthInward = location.LeakageLengthInwardsPl4;
                }
                headPlDikeTopRiver = (location.WaterLevelRiver - waterLevelRiverAverage) / (1 + (LeakageLengthOutward / LeakageLengthInward)) + pl2level;
                double deltaHead = headPlDikeTopRiver - pl2level;
                plLine.Points.Add(new GeometryPoint(firstX, 0, headValue));
                plLine.Points.Add(new GeometryPoint(dikeToeAtRiver.X, 0, headValue));
                plLine.Points.Add(new GeometryPoint(dikeTopAtRiver.X, 0, headPlDikeTopRiver));

                // calculate value at diketop polder
                double deltaXPl2 = dikeTopAtPolder.X - firstX;
                double pl2LevelDikeTopPolder = headInPlLine2Outwards + pl2gradient * deltaXPl2;
                double deltaX = dikeTopAtPolder.X - dikeTopAtRiver.X;
                double headPlDikeTopPolder = (deltaHead) * Math.Pow(Math.E, (-deltaX / LeakageLengthInward)) + pl2LevelDikeTopPolder;
                plLine.Points.Add(new GeometryPoint(dikeTopAtPolder.X, 0, headPlDikeTopPolder));

                // calculate value at dike toe
                deltaXPl2 = dikeToeAtPolder.X - firstX;
                double pl2LevelDikeToePolder = headInPlLine2Outwards + pl2gradient * deltaXPl2;
                deltaX = dikeToeAtPolder.X - dikeTopAtRiver.X;
                double headPlDikeToePolder = (deltaHead) * Math.Pow(Math.E, (-deltaX / LeakageLengthInward)) + pl2LevelDikeToePolder;
                plLine.Points.Add(new GeometryPoint(dikeToeAtPolder.X, 0, headPlDikeToePolder));

                // calculate value at 20 m after dike toe (this is needed if no uplift occurs)
                double toeplus20 = dikeToeAtPolder.X + 20;
                deltaXPl2 = toeplus20 - firstX;
                pl2LevelDikeToePolder = headInPlLine2Outwards + pl2gradient * deltaXPl2;
                deltaX = toeplus20 - dikeTopAtRiver.X;
                HeadAtDikeToePolderPlus20M = (deltaHead) * Math.Pow(Math.E, (-deltaX / LeakageLengthInward)) + pl2LevelDikeToePolder;

                // calculate value at the end of the geometry (this is needed if no uplift occurs)
                deltaXPl2 = lastX - firstX;
                pl2LevelDikeToePolder = headInPlLine2Outwards + pl2gradient * deltaXPl2;
                deltaX = lastX - dikeTopAtRiver.X;
                HeadAtEndOfGeometry = (deltaHead) * Math.Pow(Math.E, (-deltaX / LeakageLengthInward)) + pl2LevelDikeToePolder;

                // now calculate tail from slope with formula
                AddtailOfPnLine(firstX, dikeTopAtRiver.X, location, plLine, deltaHead, headInPlLine2Outwards, pl2gradient,
                                LeakageLengthInward);

                // Now continue PLline to the end with a slope of slopegradient
                // AddTailOfPl3OrPl4WithSlopeGradient(location, plLine);

                if (location.AdjustPl3And4ForUplift)
                {
                    var plLineWithoutNaN = new GeometryPointString();
                    var pnts = plLine.Points.Where(p => !double.IsNaN(p.X) && !double.IsNaN(p.Z));
                    if (pnts.Any())
                    {
                        plLineWithoutNaN.Points.AddRange(pnts);
                        AdjustLineDueToUplift(plLineWithoutNaN, plLineType, location);
                        plLineWithoutNaN.CondensePoints();
                        return plLineWithoutNaN;
                    }
                }

                //  EnsureDescendingLine(plLine);
                //  RemoveRedundantPoints(plLine);
            }

            plLine.CondensePoints();

            return plLine;
        }

        /// <summary>
        /// Create PL3 or PL4 for Ringtoets WTI 2017, according to the memo of A. van Duinen, for "Clay dike on clay" and "Sand dike on clay" (cases 1A and 2A). 
        /// For more info, see:  https://repos.deltares.nl/repos/delftgeosystems/delftgeosystems/trunk/documents/Geotechnics
        /// </summary>
        /// <param name="headOutside"></param>
        /// <param name="plLineType"></param>
        /// <param name="location"></param>
        /// <returns></returns>
        private GeometryPointString CreatePLLine3Or4ForRingtoetsWti2017ForDikesOnClay(double headOutside, PLLineType plLineType, Location location)
        {
            var plLine = new GeometryPointString();
            // short-cuts
            double firstX = location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X;
            double lastX = location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X;
            double headInPlLine2Outwards = location.CorrectedHeadInPLLine2Outwards;
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeToeAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
            GeometryPoint dikeTopAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint dikeToeAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);

            // Soil profile is used to check if there is really an aquifer below toe of dike at river.
            // The assumption is made that if there is an aquifer at the riverside and at the polderside,
            // that there is a connection between these aquifers.
            // In the uplift calculation there will also be a check on the existence of an aquifer. 
            if (!location.Surfaceline.HasAnnotation(CharacteristicPointType.DikeToeAtRiver))
            {
                throw new System.Exception("Characteristic point \"dike toe at river\" is not defined.");
            }
            if (!location.Surfaceline.HasAnnotation(CharacteristicPointType.DikeTopAtRiver))
            {
                throw new System.Exception("Characteristic point \"dike top at river\" is not defined.");
            }
            if (!location.Surfaceline.HasAnnotation(CharacteristicPointType.DikeTopAtPolder))
            {
                throw new System.Exception("Characteristic point \"dike top at polder\" is not defined.");
            }
            if (!location.Surfaceline.HasAnnotation(CharacteristicPointType.DikeToeAtPolder))
            {
                throw new System.Exception("Characteristic point \"dike toe at polder\" is not defined.");
            }

            SoilProfile1D actualSoilProfile = location.SoilProfile2D != null
                                                  ? GetSoilProfileBelowPoint(location.SoilProfile2D, dikeTopAtPolder.X) : location.SoilProfile1D;
            if (actualSoilProfile == null)
            {
                throw new System.Exception("No soil profile is present.");
            }

            SoilLayer1D relevantAquiferLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile, location);

            if (relevantAquiferLayer != null && !double.IsNaN(location.LeakageLengthInwardsPl3))
            {
                if (double.IsNaN(headInPlLine2Outwards))
                {
                    headInPlLine2Outwards = location.WaterLevelPolder;
                }

                double waterLevelRiverAverage = !double.IsNaN(location.WaterLevelRiverAverage) ? location.WaterLevelRiverAverage : headInPlLine2Outwards;

                // First, horizontal line between outside boundary and Dike Toe At River
                plLine.Points.Add(new GeometryPoint(firstX, 0, headOutside));
                plLine.Points.Add(new GeometryPoint(dikeToeAtRiver.X, 0, headOutside));

                // Then, add head at Dike Toe At River (formule 1 in memo)
                double headPlDikeTopRiver = CalculateHeadPL3OrPL4AtTopDikeOutside(location, plLineType);
                plLine.Points.Add(new GeometryPoint(dikeTopAtRiver.X, 0, headPlDikeTopRiver));

                // Then, add head at Dike Top at Polder (formule 2)
                double headPlDikeTopPolder = CalculateHeadPL3OrPL4AtX(location, plLineType, dikeTopAtPolder.X);
                plLine.Points.Add(new GeometryPoint(dikeTopAtPolder.X, 0, headPlDikeTopPolder));

                // Then, add head at Dike Toe at Polder (formule 2)
                double headPlDikeToePolder = CalculateHeadPL3OrPL4AtX(location, plLineType, dikeToeAtPolder.X);
                plLine.Points.Add(new GeometryPoint(dikeToeAtPolder.X, 0, headPlDikeToePolder));

                // Calculate head at 20 m after dike toe (this is needed if no uplift occurs)
                double xDikeToeAtPolderPlus20M = Math.Min(dikeToeAtPolder.X + 20, lastX);
                HeadAtDikeToePolderPlus20M = CalculateHeadPL3OrPL4AtX(location, plLineType, xDikeToeAtPolderPlus20M);
                // This head will be added in procedure AdjustLineDueToUpliftForRingtoetsWti2017Method

                // Calculate head at the end of the geometry (this is needed if no uplift occurs)
                HeadAtEndOfGeometry = CalculateHeadPL3OrPL4AtX(location, plLineType, lastX);
                // This head will be added in procedure AdjustLineDueToUpliftForRingtoetsWti2017Method

                // Now calculate the tail of the PL-line, from DikeToeAtPolder, using formule 2.
                // If uplift occurs, the head is adjusted to a maximum head corresponding to an uplift potential of 1.
                AdjustLineDueToUpliftForRingtoetsWti2017Method(plLine, plLineType, location);
            }
            plLine.CondensePoints();
            return plLine;
        }

        /// <summary>
        /// Calculate head 2 at X-coordinate, according to the memo of A. van Duinen, see:
        /// https://repos.deltares.nl/repos/delftgeosystems/delftgeosystems/trunk/documents/Geotechnics
        /// </summary>
        /// <param name="location">Location class containing all waternet data's</param>
        /// <param name="x">X-coordinate of the considered point</param>
        /// <returns>Head 2 at x</returns>
        private double CalculateHeadPl2AtX(Location location, double x)
        {
            double firstX = location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X;
            double lastX = location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X;
            double pl2Gradient = (location.CorrectedHeadInPLLine2Outwards - location.CorrectedHeadInPLLine2Inwards) / (firstX - lastX);
            double pl2Level = location.CorrectedHeadInPLLine2Outwards + pl2Gradient * (x - firstX);
            return pl2Level;
        }

        /// <summary>
        /// Calculate head 3 or head 4 at X-coordinate, according to the memo of A. van Duinen (formule 2), see:
        /// https://repos.deltares.nl/repos/delftgeosystems/delftgeosystems/trunk/documents/Geotechnics
        /// </summary>
        /// <param name="location">Location class containing all waternet data's</param>
        /// <param name="plLineType">Type of PL-line: either PL3 or PL4</param>
        /// <param name="x">X-coordinate of the considered point</param>
        /// <returns>Head 3 or head 4 at X</returns>
        private double CalculateHeadPL3OrPL4AtX(Location location, PLLineType plLineType, double x)
        {
            // Determine the leakage length
            double leakageLengthInward;
            double leakageLengthOutward;
            if (plLineType == PLLineType.PL3)
            {
                leakageLengthInward = location.LeakageLengthInwardsPl3;
            }
            else
            {
                leakageLengthInward = location.LeakageLengthInwardsPl4;
            }

            // Head PL3/PL4 at Dike Top River (formule 1  in memo)
            double headPlDikeTopRiver = CalculateHeadPL3OrPL4AtTopDikeOutside(location, plLineType);

            // Head PL3/PL4 at x (formule 2  in memo)
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            double headPl2AtDikeTopAtRiver = CalculateHeadPl2AtX(location, dikeTopAtRiver.X);
            double headPl2AtX = CalculateHeadPl2AtX(location, x);
            double deltaX = x - dikeTopAtRiver.X;
            double head = (headPlDikeTopRiver - headPl2AtDikeTopAtRiver) * Math.Pow(Math.E, (-deltaX / leakageLengthInward)) + headPl2AtX;
            return head;
        }

        /// <summary>
        /// Calculate head 3 or head 4 at Dike Tope Outside, according to the memo of A. van Duinen (formule 1), see:
        /// https://repos.deltares.nl/repos/delftgeosystems/delftgeosystems/trunk/documents/Geotechnics
        /// </summary>
        /// <param name="location"></param>
        /// <param name="plLineType"></param>
        /// <returns></returns>
        private double CalculateHeadPL3OrPL4AtTopDikeOutside(Location location, PLLineType plLineType)
        {
            // Determine the leakage length
            double leakageLengthInward;
            double leakageLengthOutward;
            if (plLineType == PLLineType.PL3)
            {
                leakageLengthInward = location.LeakageLengthInwardsPl3;
                leakageLengthOutward = location.LeakageLengthOutwardsPl3;
            }
            else
            {
                leakageLengthInward = location.LeakageLengthInwardsPl4;
                leakageLengthOutward = location.LeakageLengthOutwardsPl4;
            }

            double waterLevelRiverAverage = !double.IsNaN(location.WaterLevelRiverAverage)
                                                ? location.WaterLevelRiverAverage
                                                : location.CorrectedHeadInPLLine2Outwards;

            // Head PL3/PL4 at Dike Top River (formule 1  in memo)
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            double headPl2DikeTopAtRiver = CalculateHeadPl2AtX(location, dikeTopAtRiver.X);
            double headPlDikeTopRiver = (location.WaterLevelRiver - waterLevelRiverAverage)
                                        / (1 + (leakageLengthOutward / leakageLengthInward)) + headPl2DikeTopAtRiver;
            return headPlDikeTopRiver;
        }

        /// <summary>
        /// Calculate the maximum head (of PL3 or PL4) for which no uplift will occur
        /// </summary>
        /// <param name="xCoor"></param>
        /// <param name="plLineType"></param>
        /// <param name="location"></param>
        /// <returns>maximum head for which no uplift will occur</returns>
        private double DetermineMaxHead(double xCoor, PLLineType plLineType, Location location)
        {
            double volumicWeightOfWater = Physics.UnitWeightOfwater;
            SoilProfile1D actualSoilProfile = location.SoilProfile2D != null
                                                  ? GetSoilProfileBelowPoint(location.SoilProfile2D, xCoor)
                                                  : location.SoilProfile1D;
            //   AdaptSoilProfileForSurfaceLevelAccuracy(actualSoilProfile, upliftCalculator.SurfaceLevel);
            SoilLayer1D relevantSandLayer = null;
            try
            {
                relevantSandLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile, location);
            }
            catch
            {
                return Double.MaxValue;
            }
            double topLevelAquifer = relevantSandLayer.TopLevel;
            double phreaticLevel = phreaticLine.GetZAtX(xCoor);
            // calculate total pressure (= soil mass + eventually free water above surface line)
            double totalPressure = actualSoilProfile.GetTotalPressure(topLevelAquifer, phreaticLevel, volumicWeightOfWater);
            // calculate the maximum head for which no uplift will occur (i.e. uplift potential = mass / [ (head - ztop) * Gamma_water] > cUpliftFactorEquilibrium = 1)
            double maxHead = topLevelAquifer + totalPressure / (cUpliftFactorEquilibrium * volumicWeightOfWater);
            return maxHead;
        }

        private void AdjustLineDueToUplift(GeometryPointString plLine, PLLineType plLineType, Location location)
        {
            int index = -1;
            double xstart = 0;
            double head = 0;
            double xend = 0;
            double maxHead = 0;
            double maxHeadplus20 = 0;
            double headplus20 = 0;
            double gradient = 0;
            double gradientEnd = 0;
            double UpliftLength = 20;
            double searchUpliftIncrement = 0.5;
            double startx = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
            double lastx = plLine.Points[plLine.Points.Count - 1].X;
            bool isUplift = false;
            xstart = startx - searchUpliftIncrement;

            while ((!isUplift) && ((xstart + searchUpliftIncrement) < lastx))
            {
                xstart += searchUpliftIncrement;
                maxHead = DetermineMaxHead(xstart, plLineType, location);
                head = plLine.GetZAtX(xstart);
                isUplift = (head > maxHead);
            }

            if (isUplift)
            {
                // create uplift point
                index = FindIndexOfPointInGeometryPointStringBeforeXCoordinate(plLine, xstart);
                if (index >= 0)
                {
                    if (plLine.Points[index + 1].X == xstart)
                    {
                        plLine.Points[index + 1].Z = maxHead;
                    }
                    else
                    {
                        var upliftPoint = new GeometryPoint(xstart, 0, maxHead);
                        plLine.Points.Insert(index + 1, upliftPoint);
                    }
                    xend = xstart + UpliftLength;
                    if (location.Surfaceline.HasAnnotation(CharacteristicPointType.DitchDikeSide) &&
                        location.Surfaceline.HasAnnotation(CharacteristicPointType.DitchPolderSide))
                    {
                        double xstartditch =
                            location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X;
                        double xendditch =
                            location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X;
                        // end of uplift not in ditch
                        if ((xend > xstartditch) && (xend < xendditch))
                        {
                            xend = xendditch;
                            UpliftLength = xend - xstart;
                        }
                    }
                    // uplift zone can not exceed geometry boundary
                    if (xend > lastx)
                    {
                        xend = lastx;
                        UpliftLength = xend - xstart;
                    }

                    maxHeadplus20 = DetermineMaxHead(xend, plLineType, location);
                    headplus20 = plLine.GetZAtX(xend);
                    headplus20 = Math.Min(headplus20, maxHeadplus20);
                    gradient = (headplus20 - maxHead) / UpliftLength;
                    gradientEnd = (location.CorrectedHeadInPLLine2Inwards - headplus20) / (lastx - xend);

                    index = FindIndexOfPointInGeometryPointStringBeforeXCoordinate(plLine, xend);
                    if (index >= 0)
                    {
                        if (plLine.Points[index + 1].X == xend)
                        {
                            // if the point already excists change the head 
                            plLine.Points[index + 1].Z = headplus20;
                        }
                        else
                        {
                            // if it is a new point add the point
                            var endupliftPoint = new GeometryPoint(xend, 0, headplus20);
                            plLine.Points.Insert(index + 1, endupliftPoint);
                        }
                    }
                    // bring the points in the uplift area (usual 20 m) to uplift level

                    for (int j = 0; j < plLine.Points.Count; j++)
                    {
                        if ((plLine.Points[j].X > xstart) && (plLine.Points[j].X < xend))
                        {
                            plLine.Points[j].Z = maxHead + (plLine.Points[j].X - xstart) * gradient;
                        }
                        else
                            // pl lines goe to pl2 at end geometry
                            // meaning that the calculation result depends on geometry
                            if (plLine.Points[j].X > xend)
                            {
                                plLine.Points[j].Z = headplus20 + (plLine.Points[j].X - xend) * gradientEnd;
                            }
                    }
                }
            }
            else
            {
                // if no uplift occures add point at 20 m after toe dike polder
                xstart = startx + 20;
                index = FindIndexOfPointInGeometryPointStringBeforeXCoordinate(plLine, xstart);
                if (index >= 0)
                {
                    if (!(plLine.Points[index].X == xstart))
                    {
                        var upliftPoint = new GeometryPoint(xstart, 0, HeadAtDikeToePolderPlus20M);
                        plLine.Points.Insert(index + 1, upliftPoint);
                    }

                    gradientEnd = (HeadAtEndOfGeometry - HeadAtDikeToePolderPlus20M) / (lastx - (xstart));

                    for (int k = 0; k < plLine.Points.Count; k++)
                    {
                        if (plLine.Points[k].X > (xstart))
                        {
                            plLine.Points[k].Z = HeadAtDikeToePolderPlus20M + (plLine.Points[k].X - (xstart)) * gradientEnd;
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Adjust PL3 or 4 if uplift occurs. 
        /// </summary>
        /// <param name="plLine"></param>
        /// <param name="plLineType"> </param>
        /// <param name="location"></param>
        private void AdjustLineDueToUpliftForRingtoetsWti2017Method(GeometryPointString plLine, PLLineType plLineType, Location location)
        {
            GeometryPoint ditchDikeSide = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
            GeometryPoint bottomDitchDikeSide = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchDikeSide);
            double head = 0;
            double maxHeadFirstPoint = 0;
            double xFirstPointPlus20M = 0;
            double upliftLength = 20;
            const double searchUpliftIncrement = 0.5;
            double firstX = location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X;
            double dikeToePolderX = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
            double lastX = location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X;
            bool isUplift = false;
            bool isDitchPresent = (ditchDikeSide != null && !double.IsNaN(ditchDikeSide.X)) && (bottomDitchDikeSide != null && !double.IsNaN(bottomDitchDikeSide.X));

            // First, create the list of X-coordinates used for the search area (= points of the surface line + intermediaire points every 50 cm)
            var listXCoord = new List<double>();
            for (int i = 0; i < location.SoilProfile2D.Geometry.Points.Count; i++)
            {
                if (location.SoilProfile2D.Geometry.Points[i].X > dikeToePolderX)
                {
                    listXCoord.Add(location.SoilProfile2D.Geometry.Points[i].X);
                }
            }
            for (int i = 0; i < (lastX - dikeToePolderX) / searchUpliftIncrement - 1; i++)
            {
                listXCoord.Add(dikeToePolderX + (i + 1) * searchUpliftIncrement);
            }
            listXCoord.Sort();

            // Search first X-coordinate at which uplift occurs (only if no ditch is present)
            double xFirstPointUplift;
            xFirstPointUplift = dikeToePolderX - searchUpliftIncrement;
            while ((!isUplift) && ((xFirstPointUplift + searchUpliftIncrement) < lastX))
            {
                xFirstPointUplift += searchUpliftIncrement;
                maxHeadFirstPoint = DetermineMaxHead(xFirstPointUplift, plLineType, location);
                head = CalculateHeadPL3OrPL4AtX(location, plLineType, xFirstPointUplift);
                isUplift = (head > maxHeadFirstPoint);
            }
            // If there is uplift and a Ditch and the upliftpoint is beyond the ditch (as seen from the diketoe), then reset upliftpoint (and head) to 
            // ditch dike side.
            if (isUplift && isDitchPresent)
            {
                if (xFirstPointUplift > ditchDikeSide.X)
                {
                    xFirstPointUplift = ditchDikeSide.X;
                    maxHeadFirstPoint = DetermineMaxHead(xFirstPointUplift, plLineType, location);
                }
            }

            if (isUplift && location.AdjustPl3And4ForUplift)
            {
                // In situation with an uplift zone of 20 m, the points situated at DikeTopAtPolder and DikeToeAtPolder are deleted
                // (i.e. PL-line is a straight line between DikeTopAtRiver and the first point at which uplift occurs)
                xFirstPointPlus20M = xFirstPointUplift + upliftLength;
                if (isDitchPresent)
                {
                    xFirstPointPlus20M = ditchDikeSide.X + upliftLength;
                }
                // uplift zone can not exceed geometry boundary inside
                if (xFirstPointPlus20M > lastX)
                {
                    xFirstPointPlus20M = lastX;
                }
                double maxHeadFirstPointPlus20M = DetermineMaxHead(xFirstPointPlus20M, plLineType, location);
                double head20M = CalculateHeadPL3OrPL4AtX(location, plLineType, xFirstPointPlus20M);
                bool isUpliftAfter20M = (head20M > maxHeadFirstPointPlus20M);
                if (isUpliftAfter20M)
                {
                    plLine.Points.Remove(plLine.Points.Last());
                    plLine.Points.Remove(plLine.Points.Last());
                }
                
                // Add the first point at which uplift occurs
                double headFirstPointUplift = CalculateHeadPL3OrPL4AtX(location, plLineType, xFirstPointUplift);
                headFirstPointUplift = Math.Min(headFirstPointUplift, maxHeadFirstPoint);

                // When the first location where uplift occurs coincides with Dike Toe at Polder, this point replaces the current point of the PL
                if (Math.Abs(xFirstPointUplift - dikeToePolderX) < 0.0001)
                {
                    plLine.Points.Remove(plLine.Points.Last());
                }
                plLine.Points.Add(new GeometryPoint(xFirstPointUplift, 0, headFirstPointUplift));

                // Add the first point at which uplift occurs + 20 m 
                double headFirstPointPlus20M = CalculateHeadPL3OrPL4AtX(location, plLineType, xFirstPointPlus20M);
                headFirstPointPlus20M = Math.Min(headFirstPointPlus20M, maxHeadFirstPointPlus20M);
                plLine.Points.Add(new GeometryPoint(xFirstPointPlus20M, 0, headFirstPointPlus20M));

                // Fill list with all relevant points in uplift area
                List<double> xValues = new List<double>();
                IList<GeometryPoint> points = location.SoilProfile2D != null ? location.SoilProfile2D.Geometry.Points : location.Surfaceline.Geometry.Points;
                foreach (var point in points)
                {
                    if (!double.IsNaN(point.X) && point.X > xFirstPointUplift && point.X < xFirstPointPlus20M && !xValues.Contains(Math.Round(point.X, 3)))
                    {
                        xValues.Add(Math.Round(point.X, 3));
                    }
                }
                xValues.Add(xFirstPointUplift);
                xValues.Add(xFirstPointPlus20M);
                xValues.Sort();

                // Fill large gaps between points
                points.Sort();
                int pointIndex = 0;
                while (pointIndex < xValues.Count - 1)
                {
                    if (xValues[pointIndex + 1] - xValues[pointIndex] > searchUpliftIncrement)
                    {
                        xValues.Insert(pointIndex + 1, Math.Min(xValues[pointIndex] + searchUpliftIncrement, (xValues[pointIndex] + xValues[pointIndex + 1]) / 2));
                    }

                    pointIndex++;
                }
                
                // Add all the points of the surface line part of the uplift zone (i.e. from first point at which uplift occurs until 20 m)
                foreach (double x in xValues)
                {
                    double xCoor = x;
                    double headAtX = CalculateHeadPL3OrPL4AtX(location, plLineType, xCoor);
                    double maxHeadAtX = DetermineMaxHead(xCoor, plLineType, location);
                    headAtX = Math.Min(headAtX, maxHeadAtX);
                    plLine.Points.Add(new GeometryPoint(xCoor, 0, headAtX));
                }
                
                // Now continue PL-line to PL2 (or PL1 if PL2 is not present), with a slope of 1:100.
                plLine.Points.Sort();
                plLine.RemoveUnnecessaryPoints();

                GeometryPoint lastCurrentPointPlLine = plLine.Points.Last();
                const double slopeEnd = 100;
                if (Math.Abs(plLine.Points.Last().X - lastX) > 0.0001)
                {
                    double head2Or1AtLastX;
                    if (Pl2 != null)
                    {
                        head2Or1AtLastX = CalculateHeadPl2AtX(location, lastX);
                    }
                    else
                    {
                        head2Or1AtLastX = phreaticLine.GetZAtX(lastX);
                    }
                    GeometryPoint pointE = DetermineIntersectionBetweenPl2Or1AndFinalSlope100OfPl3(lastCurrentPointPlLine, slopeEnd, location);
                    if (pointE.X > lastX)
                    {
                        // 1 point is added (on the slope 1:100) at end of the geometry
                        double headOnSlopeAtLastX;
                        bool slopeGoesDown = (CalculateHeadPl2AtX(location, lastCurrentPointPlLine.X) < lastCurrentPointPlLine.Z);
                        if (slopeGoesDown)
                        {
                            headOnSlopeAtLastX = lastCurrentPointPlLine.Z - (lastX - lastCurrentPointPlLine.X) / 100;
                        }
                        else
                        {
                            headOnSlopeAtLastX = lastCurrentPointPlLine.Z + (lastX - lastCurrentPointPlLine.X) / 100;
                        }                     
                        plLine.Points.Add(new GeometryPoint(lastX, 0, headOnSlopeAtLastX));
                    }
                    else
                    {
                        // 2 points are added: 1 at the end of the slope 1:100 (point E in figure 7 of the memo of A. van Duinen) and 1 at end of geometry (point G)
                        if (pointE != null)
                        {
                            plLine.Points.Add(pointE);
                        }
                        plLine.Points.Add(new GeometryPoint(lastX, 0, head2Or1AtLastX));
                    }
                }
            }
            else
            {
                // if no uplift occurs, add point at 20 m after dike toe at polder
                xFirstPointPlus20M = Math.Min(dikeToePolderX + 20, lastX);
                plLine.Points.Add(new GeometryPoint(xFirstPointPlus20M, 0, HeadAtDikeToePolderPlus20M));
                // and add point at the end of geometry
                plLine.Points.Add(new GeometryPoint(lastX, 0, HeadAtEndOfGeometry));
            }
            plLine.Points.Sort();
        }

        /// <summary>
        ///     The correction is applied as follows:
        ///     Check every point from DikeToeAtPolder to SurfaceLevelInside (from left to right) for uplift.
        ///     If uplift occurs, then correct PL3/PL4 value, so uplift does not occur.
        ///     All points in PL3 from this point to DikeToeAtRiver should be removed.
        ///     The PL3 continues from this point on with the specified slopegradient until polderlevel.
        ///     Make sure PL3 is always descending from left to right.
        ///     A better implementation (not implemented yet) would be:
        ///     Correct plline 3 or 4 for uplift according to
        ///     TRW (Technisch Rapport Waterspanningen bij dijken) par. b1.3.4 "Stijghoogte in het eerste watervoerende pakket"
        ///     - Adjust PL3/4 for all surface points from end of profile to toe of dike, so no uplift will occur in that surface
        ///     point
        ///     - From the point, closest to the dike, (firstAdjustedPLPoint) where this correction has been made the following has
        ///     to be done
        ///     * PL3/4 will continue horizontally from firstAdjustedPLPoint over a distance L = 2* d (d is height all layers above
        ///     the aquifer)
        ///     * The the PL3/4 will go down in a slope of 1:50 to the PolderLevel
        ///     PL3/4-----
        ///               \___________ L = 2 * d
        ///                           \
        ///                            \__________
        /// </summary>
        /// <param name="plLine"></param>
        /// <param name="plLineType"> </param>
        /// <param name="location"></param>
        private void AdjustLineAccordingToTRWUplift(GeometryPointString plLine, PLLineType plLineType, Location location)
        {
            ClearOutputValuesForPl3_4(plLineType);
            var upliftCalculator = new UpliftCalculator();

            upliftCalculator.IsUseOvenDryUnitWeight = IsUseOvenDryUnitWeight;
            GeometryPoint startSurfacePoint = GetDikeToeInward(location.Surfaceline);
            int indexOfFixedPlPoint = FindIndexOfPointInGeometryPointStringWithXCoordinate(plLine,
                                                                                           location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X);
            if (indexOfFixedPlPoint < 0)
            {
                throw new System.Exception("Could not find fixed point for PLLine");
            }

            // Adjust PL3/4 for all surface points from toe of dike to end of profile to, so no uplift will occur in that surface point
            IEnumerable<GeometryPoint> relevantSurfacePointsList =
                location.Surfaceline.Geometry.Points.Cast<GeometryPoint>()
                        .Where(point => point.X >= startSurfacePoint.X)
                        .OrderBy(point => point.X);
            // Adjustment will only be applied if the value to adjust to is smaller than the previous adjusted valuet (to avoid that the PL3/PL4 will be adjusted to the 
            // polderside of the ditch i.s.o. the dikeside of the ditch.
            // So we remember which was the last adjsuted value in lastAdjustedHeadOfPlLine.
            double lastAdjustedHeadOfPlLine = Double.MaxValue;
            foreach (var surfacePoint in relevantSurfacePointsList)
            {
                upliftCalculator.SurfaceLevel = surfacePoint.Z;
                if (phreaticLine != null)
                {
                    upliftCalculator.PhreaticLevel = LineHelper.ZFromX(surfacePoint.X, phreaticLine.Points);
                    // set phreatic level to calculate upliftfactor
                }
                else
                {
                    upliftCalculator.PhreaticLevel = location.WaterLevelPolder;
                    // if not PL1 created then assume phreatic level is same as waterlevel polder
                }
                double headOfPLLine = LineHelper.ZFromX(surfacePoint.X, plLine.Points);
                SoilProfile1D actualSoilProfile = location.SoilProfile2D != null
                                                      ? GetSoilProfileBelowPoint(location.SoilProfile2D, surfacePoint.X)
                                                      : location.SoilProfile1D;
                AdaptSoilProfileForSurfaceLevelAccuracy(actualSoilProfile, upliftCalculator.SurfaceLevel);
                SoilLayer1D relevantSandLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile, location);
                // When plLineType is PL4 it is possible that no relevant aquifer layer is found, then the adjustment is not necessary 
                if (relevantSandLayer != null)
                {
                    upliftCalculator.TopOfLayerToBeEvaluated = relevantSandLayer.TopLevel;
                    upliftCalculator.SoilProfile = actualSoilProfile;
                    double upliftFactor = upliftCalculator.CalculateUpliftFactor(headOfPLLine);
                    if (upliftFactor < cUpliftFactorEquilibrium)
                    {
                        // Adjust headOfPLLine so there is equilibrium
                        headOfPLLine = Math.Max(upliftCalculator.CalculateHeadOfPLLine(cUpliftFactorEquilibrium),
                                                waterLevelPolder);
                        //adjustments according to changes in DAM
                        if (headOfPLLine < lastAdjustedHeadOfPlLine)
                        {
                            var plPoint = new GeometryPoint(surfacePoint.X, 0, headOfPLLine);
                            // Remove all points of PLLine after fixed point
                            LineHelper.RemoveAllPointsOfGeometryPointStringAfterIndex(plLine, indexOfFixedPlPoint);
                            plLine.Points.Add(plPoint);
                            AddTailOfPl3OrPl4WithSlopeGradient(location, plLine);
                            lastAdjustedHeadOfPlLine = headOfPLLine;
                            UpdateOutputValuesForPl3_4(plLineType, surfacePoint.X, upliftFactor, headOfPLLine);
                        }
                    }
                    UpdateOutputValuesForPl3_4(plLineType, surfacePoint.X, upliftFactor, headOfPLLine);
                }
            }
        }

        private void ClearOutputValuesForPl3_4(PLLineType plLineType)
        {
            switch (plLineType)
            {
                case PLLineType.PL3:
                    pl3MinUplift = Double.MaxValue;
                    pl3HeadAdjusted = Double.MaxValue;
                    pl3LocationXMinUplift = Double.MaxValue;
                    upliftValuesPl3 = new Dictionary<double, double>();
                    break;
                case PLLineType.PL4:
                    pl4MinUplift = Double.MaxValue;
                    pl4HeadAdjusted = Double.MaxValue;
                    pl4LocationXMinUplift = Double.MaxValue;
                    upliftValuesPl4 = new Dictionary<double, double>();
                    break;
            }
        }

        private void UpdateOutputValuesForPl3_4(PLLineType plLineType, double xCoordinate, double upliftFactor,
                                                double headValue)
        {
            if (Math.Abs(upliftFactor - double.MaxValue) > 100)
            {
                switch (plLineType)
                {
                    case PLLineType.PL3:
                        pl3MinUplift = upliftFactor;
                        pl3HeadAdjusted = headValue;
                        pl3LocationXMinUplift = xCoordinate;
                        if (!upliftValuesPl3.ContainsKey(xCoordinate))
                        {
                            upliftValuesPl3[xCoordinate] = upliftFactor;
                        }
                        break;
                    case PLLineType.PL4:

                        pl4MinUplift = upliftFactor;
                        pl4HeadAdjusted = headValue;
                        pl4LocationXMinUplift = xCoordinate;
                        upliftValuesPl4[xCoordinate] = upliftFactor;
                        break;
                }
            }
        }

        /// <summary>
        ///     Finds the index of point with X coordinate.
        /// </summary>
        /// <param name="plLine">The pl line.</param>
        /// <param name="x">The x coordinate.</param>
        /// <returns></returns>
        private int FindIndexOfPointInGeometryPointStringWithXCoordinate(GeometryPointString plLine, double x)
        {
            for (int pointIndex = plLine.Points.Count - 1; pointIndex > 0; pointIndex--)
            {
                GeometryPoint plPoint = plLine.Points[pointIndex];
                if (plPoint.X.AlmostEquals(x, GeometryPoint.Precision))
                {
                    return pointIndex;
                }
            }

            return -1;
        }

        /// <summary>
        /// Find the index of the point before the z coordinate. If no point is found return -1
        /// </summary>
        /// <param name="plLine"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        private int FindIndexOfPointInGeometryPointStringBeforeXCoordinate(GeometryPointString plLine, double x)
        {
            for (int pointIndex = plLine.Points.Count - 1; pointIndex > 0; pointIndex--)
            {
                GeometryPoint plPoint = plLine.Points[pointIndex];
                if (plPoint.X < x)
                {
                    return pointIndex;
                }
            }
            return -1;
        }

        /// <summary>
        /// Due to accuracy surfacelevel could be slightly below toplevel of toplevel soilprofile
        /// Adjust toplevel of soilprofile to avoid this
        /// </summary>
        /// <param name="actualSoilProfile"></param>
        /// <param name="surfaceLevel"></param>
        private static void AdaptSoilProfileForSurfaceLevelAccuracy(SoilProfile1D actualSoilProfile, double surfaceLevel)
        {
            if (surfaceLevel > actualSoilProfile.TopLevel)
            {
                actualSoilProfile.Layers[0].TopLevel = surfaceLevel;
            }
        }

        private void AddtailOfPnLine(double firstX, double xTopRiver, Location location,
                                     GeometryPointString plLine, double deltaHead, double startheadPl2, double pl2gradient, double LeakageLength)
        // calculate plvalue at each geometry point
        {
            double xstart = plLine.Points[plLine.Points.Count - 1].X;
            for (int i = 0; i < location.Surfaceline.Geometry.Count; i++)
            {
                if ((!double.IsNaN(location.Surfaceline.Geometry.Points[i].X)) && (location.Surfaceline.Geometry.Points[i].X > xstart))
                {
                    double xcoor = location.Surfaceline.Geometry.Points[i].X;
                    double deltaXPl2 = xcoor - firstX;
                    double Pl2x = startheadPl2 + pl2gradient * deltaXPl2;

                    double deltaX = xcoor - xTopRiver;
                    double HeadPl = (deltaHead) * Math.Pow(Math.E, (-deltaX / LeakageLength)) + Pl2x;
                    plLine.Points.Add(new GeometryPoint(location.Surfaceline.Geometry.Points[i].X, 0, HeadPl));
                }
            }
        }

        /// <summary>
        ///     Continue PLline to the end with a slope of slopegradient
        ///     If PLLine is lower then polderlevel, continue with polderlevel
        /// </summary>
        /// <param name="location"></param>
        /// <param name="plLine"></param>
        private void AddTailOfPl3OrPl4WithSlopeGradient(Location location, GeometryPointString plLine)
        {
            if (plLine.Points.Last().Z <= location.WaterLevelPolder)
            {
                // the PL line is already at WaterLevelPolder, so countinue it horizontally
                plLine.Points.Add(new GeometryPoint(location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X, 0,
                                                    location.WaterLevelPolder));
            }
            else
            {
                //location.Surfaceline.Points.Last(p => !double.IsNaN(p.X)).X
                //var lastPoint = location.Surfaceline.Points.Max(p => p.X);
                double lastPointX = location.Surfaceline.Geometry.Points.Max(p => p.X);
                double lengthFromLastPlPointToEnd = Math.Abs(lastPointX - plLine.Points.Last(p => !double.IsNaN(p.X)).X);
                double headAtLastPlPoint = plLine.Points.Last(p => !double.IsNaN(p.Z)).Z;
                double headAtEnd = headAtLastPlPoint -
                                   location.SlopeDampingPiezometricHeightPolderSide * lengthFromLastPlPointToEnd;
                Line waterLevelPolderLine =
                    CreateNewLine(
                        new GeometryPoint(location.Surfaceline.Geometry.Points.First(p => !double.IsNaN(p.X)).X, 0,
                                          location.WaterLevelPolder),
                        new GeometryPoint(location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X, 0,
                                          location.WaterLevelPolder));
                Line slopeLine = CreateNewLine(new GeometryPoint(
                                                   location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, 0,
                                                   headAtLastPlPoint),
                                               new GeometryPoint(lastPointX, 0, headAtEnd));
                GeometryPoint intersectionPoint = LineHelper.GetIntersectPoint(waterLevelPolderLine, slopeLine,
                                                                               PointType.XZ);
                if (!double.IsNaN(intersectionPoint.X))
                {
                    plLine.Points.Add(new GeometryPoint(intersectionPoint.X, 0, location.WaterLevelPolder));
                    plLine.Points.Add(new GeometryPoint(lastPointX, 0, location.WaterLevelPolder));
                }
                else
                {
                    plLine.Points.Add(new GeometryPoint(location.Surfaceline.Geometry.Points.Last(p => !double.IsNaN(p.X)).X, 0,
                                                        headAtEnd));
                }
            }
        }

        private static SoilLayer1D GetRelevantAquiferLayer(PLLineType type, SoilProfile1D actualSoilProfile, Location location)
        {
            IList<SoilLayer1D> aquiferLayers = actualSoilProfile.GetAquiferLayers();

            SoilLayer1D relevantAquiferLayer = null;
            switch (type)
            {
                case PLLineType.PL3:
                    if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017 && location.DikeSoilScenario == DikeSoilScenario.ClayDikeOnSand)
                    {
                        relevantAquiferLayer = actualSoilProfile.TopAquiferLayer; // Highest aquifer
                    }
                    else
                    {
                        relevantAquiferLayer = actualSoilProfile.BottomAquiferLayer; // Pleistocene
                    }
                    break;
                case PLLineType.PL4:
                    relevantAquiferLayer = actualSoilProfile.InBetweenAquiferLayer; // Intermediate sand layer 
                    break;
                default:
                    string message = String.Format("Invalid PL line type:{0} for creation of PL Line 3 or 4", type);
                    throw new System.Exception(message);
            }

            if (relevantAquiferLayer == null)
            {
                string message = "Soil profile (" + actualSoilProfile.Name + ") contains no relevant aquifer layer.";
                throw new System.Exception(message);
            }
            return relevantAquiferLayer;
        }

        /// <summary>
        /// Remove redundant points (i.e. lying on a line between its both neighbours)
        /// </summary>
        /// <param name="Tolerance"></param>
        /// <param name="plLine"></param>
        private static void RemoveRedundantPoints(GeometryPointString plLine)
        {
            const double tolerance = 0.001;
            for (int pointIndex = plLine.Points.Count - 2; pointIndex > 0; pointIndex--)
            {
                GeometryPoint prev = plLine.Points[pointIndex - 1];
                GeometryPoint point = plLine.Points[pointIndex];
                GeometryPoint next = plLine.Points[pointIndex + 1];
                if (Math.Abs((point.Z - prev.Z) / (point.X - prev.X) - (next.Z - point.Z) / (next.X - point.X)) < tolerance)
                {
                    plLine.Points.Remove(point);
                }
            }
        }

        /// <summary>
        ///     All points should have descending z from dike to polder
        /// </summary>
        /// <param name="plLine"></param>
        private static void EnsureDescendingLine(GeometryPointString plLine)
        {
            GeometryPoint plPointPrevious = null;
            foreach (var point in plLine.Points)
            {
                if (plPointPrevious != null && point.Z > plPointPrevious.Z)
                {
                    point.Z = plPointPrevious.Z;
                }
                plPointPrevious = point;
            }
        }

        private IEnumerable<GeometryPoint> FindAllPlLinePointsAtNWOLocation(GeometryPointString pl1Line,
                                                                            SurfaceLine2 surfaceLine)
        {
            IEnumerable<GeometryPoint> results =
                LineHelper.GetPointSegmentBetween(
                    surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X -
                    cOffsetPhreaticLineBelowSurface,
                    surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).X +
                    cOffsetPhreaticLineBelowSurface, pl1Line.Points);

            return results;
        }

        private void AdaptPL1ForNonWaterRetainingObject(GeometryPointString pl1Line, Location location)
        {
            // check if nwo points are available as CharacteristicPoints
            GeometryPoint nwo1 = null;
            if (location.Surfaceline.IsDefined(CharacteristicPointType.NonWaterRetainingObjectPoint1))
            {
                nwo1 = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1);
            }
            GeometryPoint nwo2 = null;
            if (location.Surfaceline.IsDefined(CharacteristicPointType.NonWaterRetainingObjectPoint2))
            {
                nwo2 = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint2);
            }
            GeometryPoint nwo3 = null;
            if (location.Surfaceline.IsDefined(CharacteristicPointType.NonWaterRetainingObjectPoint3))
            {
                nwo3 = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint3);
            }
            GeometryPoint nwo4 = null;
            if (location.Surfaceline.IsDefined(CharacteristicPointType.NonWaterRetainingObjectPoint4))
            {
                nwo4 = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4);
            }
            if ((nwo1 != null) && (nwo2 != null) && (nwo3 != null) && (nwo4 != null))
            {
                // Find all points in the Pl line that (might) coincide with the NWO
                IEnumerable<GeometryPoint> plPointsToBeMoved = FindAllPlLinePointsAtNWOLocation(pl1Line,
                                                                                                location.Surfaceline);

                GeometryPoint nwo1Pl = null;
                GeometryPoint nwo2Pl = null;
                GeometryPoint nwo3Pl = null;
                GeometryPoint nwo4Pl = null;

                // For NWO point, determine whether a pl point has to be added 
                if (LineHelper.PositionXzOfPointRelatedToPLLine(nwo1, pl1Line.Points) !=
                    PLLinePointPositionXzType.AbovePLLine)
                {
                    nwo1Pl = new GeometryPoint(nwo1.X, 0, nwo1.Z - cOffsetPhreaticLineBelowSurface);
                }
                if (LineHelper.PositionXzOfPointRelatedToPLLine(nwo2, pl1Line.Points) !=
                    PLLinePointPositionXzType.AbovePLLine)
                {
                    nwo2Pl = new GeometryPoint(nwo2.X, 0, nwo2.Z - cOffsetPhreaticLineBelowSurface);
                }
                if (LineHelper.PositionXzOfPointRelatedToPLLine(nwo3, pl1Line.Points) !=
                    PLLinePointPositionXzType.AbovePLLine)
                {
                    nwo3Pl = new GeometryPoint(nwo3.X, 0, nwo3.Z - cOffsetPhreaticLineBelowSurface);
                }
                if (LineHelper.PositionXzOfPointRelatedToPLLine(nwo4, pl1Line.Points) !=
                    PLLinePointPositionXzType.AbovePLLine)
                {
                    nwo4Pl = new GeometryPoint(nwo4.X, 0, nwo4.Z - cOffsetPhreaticLineBelowSurface);
                }

                // Find the intersections of pl line and NWO and handle these
                // Intersection between nwo point1 and nwo point2 only when nwo point1 is above pl line and nwo point2 is below plLine
                GeometryPoint intersection1 = null;
                if ((nwo1Pl == null) && (nwo2Pl != null))
                {
                    Line lineNWO = CreateNewLine(nwo1, nwo2);
                    IList<GeometryPoint> ips = pl1Line.IntersectionPointsXzWithLineXz(lineNWO);
                    if (ips.Count > 0)
                    {
                        intersection1 = new GeometryPoint(ips.First().X, 0, ips.First().Z);
                    }
                }
                // Intersection between nwo point3 and nwo point4 only when nwo point3 is below pl line and nwo point4 is above plLine
                GeometryPoint intersection2 = null;
                if ((nwo3Pl != null) && (nwo4Pl == null))
                {
                    Line lineNWO = CreateNewLine(nwo3, nwo4);
                    IList<GeometryPoint> ips = pl1Line.IntersectionPointsXzWithLineXz(lineNWO);
                    if (ips.Count > 0)
                    {
                        intersection2 = new GeometryPoint(ips.Last().X, 0, ips.Last().Z);
                    }
                }

                // Handle making the NWO empty
                if ((location.NWOPhreaticAdaption != PhreaticAdaptionType.None) &&
                    (location.NWOPhreaticAdaption == PhreaticAdaptionType.MakeEmpty))
                {
                    // for the polderside, the pl line is always allowed to be adapted. For the riverside, the pl line may only be adapted when the original waterlevel is runs through the nwo.
                    RemoveAllWaterFromNonWaterRetainingObject(nwo1, pl1Line, nwo1Pl, nwo2Pl, nwo3Pl, nwo4Pl,
                                                              intersection1, intersection2, plPointsToBeMoved,
                                                              location.Surfaceline);
                }
                // Handle making the waterlevel horizontal in the NWO at the Riverside when needed (Polderside is already done when needed, see CreatePhreaticLineSegmentsInShoulderAndPolder.   
                if (location.NWOPhreaticAdaption != PhreaticAdaptionType.None)
                {
                    // For the riverside, the pl line may only be adapted when the original waterlevel is runs through the nwo and is not already level.
                    if ((nwo1.X <= location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X) &&
                        ((intersection1 != null) || (intersection2 != null)))
                    {
                        double requiredWaterLevel;
                        // Check whether adaption of intersection points is needed
                        if (intersection2 == null)
                        {
                            // only intersection 1 avaialable so add intersection 2
                            // first see if nwo3/4 intersects, if not try nwo2/3. If still no intersection found valid level not possible, raise error
                            MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1(nwo2, nwo3, nwo4, pl1Line, nwo3Pl,
                                                                                       nwo4Pl, intersection1);
                            requiredWaterLevel = intersection1.Z;
                        }
                        else
                        {
                            if (intersection1 == null)
                            {
                                // only intersection 2 avaialable so add intersection 1
                                // first see if nwo1/2 intersects, if not try nwo2/3. If still no intersection found valid level not possible, raise error
                                MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2(nwo1, nwo2, nwo3, pl1Line,
                                                                                           nwo1Pl, nwo2Pl, nwo3Pl,
                                                                                           intersection2);
                                requiredWaterLevel = intersection2.Z;
                            }
                            else
                            {
                                // intersection 1 and intersection 2 available. Only act when levels were different.
                                requiredWaterLevel = Math.Min(intersection1.Z, intersection2.Z);
                                if ((Math.Abs(intersection1.Z - intersection2.Z) > GeometryPoint.Precision))
                                {
                                    if (intersection1.Z < intersection2.Z)
                                    {
                                        // make level in NWO intersection1.Z
                                        MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1And2(nwo2, nwo3, nwo4,
                                                                                                       pl1Line, nwo3Pl,
                                                                                                       intersection1,
                                                                                                       intersection2);
                                    }
                                    else
                                    {
                                        // make level in NWO intersection2.Z
                                        MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2And1(nwo1, nwo2, nwo3,
                                                                                                       pl1Line, nwo2Pl,
                                                                                                       intersection1,
                                                                                                       intersection2);
                                    }
                                }
                            }
                        }

                        // Move all the points in the pl line itself that need to be moved to the horizontal proper level.
                        foreach (var plLinePoint in plPointsToBeMoved)
                        {
                            plLinePoint.Z = requiredWaterLevel;
                        }
                        LineHelper.DeleteCoinsidingPoints(pl1Line.Points, GeometryPoint.Precision);
                    }
                }
            }
        }

        private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2And1(GeometryPoint nwo1,
                                                                                    GeometryPoint nwo2,
                                                                                    GeometryPoint nwo3,
                                                                                    GeometryPointString pl1Line,
                                                                                    GeometryPoint nwo2Pl,
                                                                                    GeometryPoint intersection1,
                                                                                    GeometryPoint intersection2)
        {
            Line lineNWO = CreateNewLine(nwo1, nwo2);
            Line linePL = CreateNewLine(new GeometryPoint(nwo1.X, 0, intersection2.Z), intersection2);
            GeometryPoint isp = LineHelper.GetIntersectPoint(lineNWO, linePL, PointType.XZ);
            if (!double.IsNaN(isp.X))
            {
                CreateNewPointOnPlLine(pl1Line, intersection2);
                CreateNewPointOnPlLine(pl1Line, isp, -cOffsetPhreaticLineBelowSurface, intersection2);
                CreateNewPointOnPlLine(pl1Line, intersection1, -cOffsetPhreaticLineBelowSurface);
            }
            else
            {
                Line lineNWOb = CreateNewLine(nwo2, nwo3);
                isp = LineHelper.GetIntersectPoint(lineNWOb, linePL, PointType.XZ);
                if (!double.IsNaN(isp.X))
                {
                    CreateNewPointOnPlLine(pl1Line, intersection2);
                    CreateNewPointOnPlLine(pl1Line, isp, -cOffsetPhreaticLineBelowSurface, intersection2);
                    CreateNewPointOnPlLine(pl1Line, nwo2Pl);
                    if (nwo2Pl.X > intersection1.X)
                    {
                        CreateNewPointOnPlLine(pl1Line, intersection1);
                    }
                }
                else
                {
                    string message =
                        "Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.";
                    throw new System.Exception(message);
                }
            }
        }

        private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1And2(GeometryPoint nwo2,
                                                                                    GeometryPoint nwo3,
                                                                                    GeometryPoint nwo4,
                                                                                    GeometryPointString pl1Line,
                                                                                    GeometryPoint nwo3Pl,
                                                                                    GeometryPoint intersection1,
                                                                                    GeometryPoint intersection2)
        {
            var lineNWO = new Line
                {
                    BeginPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z),
                    EndPoint = new GeometryPoint(nwo4.X, 0, nwo4.Z)
                };
            var linePL = new Line
                {
                    BeginPoint = new GeometryPoint(intersection1.X, 0, intersection1.Z),
                    EndPoint = new GeometryPoint(nwo4.X, 0, intersection1.Z)
                };
            GeometryPoint isp = LineHelper.GetIntersectPoint(lineNWO, linePL, PointType.XZ);
            if (!double.IsNaN(isp.X))
            {
                GeometryPoint newP1 = LineHelper.EnsurePointAtX(intersection1.X, pl1Line.Points, GeometryPoint.Precision);
                newP1.Z = intersection1.Z;
                GeometryPoint newP2 = LineHelper.EnsurePointAtX(isp.X + cOffsetPhreaticLineBelowSurface, pl1Line.Points,
                                                                GeometryPoint.Precision);
                newP2.Z = intersection1.Z;
                GeometryPoint newP3 = LineHelper.EnsurePointAtX(intersection2.X + cOffsetPhreaticLineBelowSurface,
                                                                pl1Line.Points,
                                                                GeometryPoint.Precision);
                newP3.Z = intersection2.Z;
            }
            else
            {
                var lineNWOb = new Line
                    {
                        BeginPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z),
                        EndPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z)
                    };

                isp = LineHelper.GetIntersectPoint(lineNWOb, linePL, PointType.XZ);
                if (!double.IsNaN(isp.X))
                {
                    GeometryPoint newP1 = LineHelper.EnsurePointAtX(intersection1.X, pl1Line.Points,
                                                                    GeometryPoint.Precision);
                    newP1.Z = intersection1.Z;
                    GeometryPoint newP2 = LineHelper.EnsurePointAtX(isp.X + cOffsetPhreaticLineBelowSurface,
                                                                    pl1Line.Points,
                                                                    GeometryPoint.Precision);
                    newP2.Z = intersection1.Z;
                    GeometryPoint newP3 = LineHelper.EnsurePointAtX(nwo3Pl.X, pl1Line.Points, GeometryPoint.Precision);
                    newP3.Z = nwo3Pl.Z;
                    if (nwo3Pl.X < intersection2.X)
                    {
                        GeometryPoint newP4 =
                            LineHelper.EnsurePointAtX(intersection2.X + cOffsetPhreaticLineBelowSurface,
                                                      pl1Line.Points, GeometryPoint.Precision);
                        newP4.Z = intersection2.Z;
                    }
                }
                else
                {
                    string message =
                        "Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.";
                    throw new System.Exception(message);
                }
            }
        }

        private void RemoveAllWaterFromNonWaterRetainingObject(GeometryPoint nwo1, GeometryPointString pl1Line,
                                                               GeometryPoint nwo1Pl, GeometryPoint nwo2Pl, GeometryPoint nwo3Pl, GeometryPoint nwo4Pl,
                                                               GeometryPoint intersection1, GeometryPoint intersection2, IEnumerable<GeometryPoint> plPointsToBeMoved,
                                                               SurfaceLine2 surfaceline)
        {
            if ((nwo1.X >= surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X) ||
                ((intersection1 != null) && (intersection2 != null)))
            {
                // Move all the points in the pl line itself that need to be moved to below the surfaceline.
                MoveSelectedPLLinePointsBelowSurfaceLine(plPointsToBeMoved, surfaceline);

                // now add all extra points to the pl line 
                if (nwo1Pl != null)
                {
                    GeometryPoint newP = LineHelper.EnsurePointAtX(nwo1Pl.X, pl1Line.Points, GeometryPoint.Precision);
                    newP.Z = nwo1Pl.Z;
                }
                if (nwo2Pl != null)
                {
                    GeometryPoint newP = LineHelper.EnsurePointAtX(nwo2Pl.X, pl1Line.Points, GeometryPoint.Precision);
                    newP.Z = nwo2Pl.Z;
                }
                if (nwo3Pl != null)
                {
                    GeometryPoint newP = LineHelper.EnsurePointAtX(nwo3Pl.X, pl1Line.Points, GeometryPoint.Precision);
                    newP.Z = nwo3Pl.Z;
                }
                if (nwo4Pl != null)
                {
                    GeometryPoint newP = LineHelper.EnsurePointAtX(nwo4Pl.X, pl1Line.Points, GeometryPoint.Precision);
                    newP.Z = nwo4Pl.Z;
                }
                // Note: for intersection points, apply offset in X direction not in Z.
                if (intersection1 != null)
                {
                    GeometryPoint newP = LineHelper.EnsurePointAtX(intersection1.X - cOffsetPhreaticLineBelowSurface,
                                                                   pl1Line.Points, GeometryPoint.Precision);
                    newP.Z = intersection1.Z;
                }
                if (intersection2 != null)
                {
                    GeometryPoint newP = LineHelper.EnsurePointAtX(intersection2.X + cOffsetPhreaticLineBelowSurface,
                                                                   pl1Line.Points, GeometryPoint.Precision);
                    newP.Z = intersection2.Z;
                }
                LineHelper.DeleteCoinsidingPoints(pl1Line.Points, GeometryPoint.Precision);
            }
        }

        private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2(GeometryPoint nwo1, GeometryPoint nwo2,
                                                                                GeometryPoint nwo3,
                                                                                GeometryPointString pl1Line,
                                                                                GeometryPoint nwo1Pl,
                                                                                GeometryPoint nwo2Pl,
                                                                                GeometryPoint nwo3Pl,
                                                                                GeometryPoint intersection2)
        {
            Line lineNWO = CreateNewLine(nwo1, nwo2);
            Line linePL = CreateNewLine(nwo1, intersection2);
            GeometryPoint isp = LineHelper.GetIntersectPoint(lineNWO, linePL, PointType.XZ);
            if (!double.IsNaN(isp.X))
            {
                CreateNewPointOnPlLine(pl1Line, intersection2);
                CreateNewPointOnPlLine(pl1Line, isp, -cOffsetPhreaticLineBelowSurface, intersection2);
                CreateNewPointOnPlLine(pl1Line, nwo1Pl);
            }
            else
            {
                Line lineNWOb = CreateNewLine(nwo2, nwo3);
                isp = LineHelper.GetIntersectPoint(lineNWOb, linePL, PointType.XZ);
                if (!double.IsNaN(isp.X))
                {
                    CreateNewPointOnPlLine(pl1Line, intersection2);
                    CreateNewPointOnPlLine(pl1Line, isp, -cOffsetPhreaticLineBelowSurface, intersection2);
                    CreateNewPointOnPlLine(pl1Line, nwo2Pl);
                    if ((nwo1Pl != null) && (nwo2Pl.X > nwo1Pl.X))
                    {
                        CreateNewPointOnPlLine(pl1Line, nwo1Pl);
                    }
                }
                else
                {
                    string message =
                        "Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.";
                    throw new System.Exception(message);
                }
            }
        }

        private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1(GeometryPoint nwo2, GeometryPoint nwo3,
                                                                                GeometryPoint nwo4,
                                                                                GeometryPointString pl1Line,
                                                                                GeometryPoint nwo3Pl,
                                                                                GeometryPoint nwo4Pl,
                                                                                GeometryPoint intersection1)
        {
            Line lineNWO = CreateNewLine(nwo3, nwo4);
            Line linePL = CreateNewLine(nwo4, intersection1);
            GeometryPoint isp = LineHelper.GetIntersectPoint(lineNWO, linePL, PointType.XZ);

            if (!double.IsNaN(isp.X))
            {
                CreateNewPointOnPlLine(pl1Line, intersection1, double.NaN, null);
                CreateNewPointOnPlLine(pl1Line, isp, cOffsetPhreaticLineBelowSurface, intersection1);
                CreateNewPointOnPlLine(pl1Line, nwo4Pl);
            }
            else
            {
                Line lineNWOb = CreateNewLine(nwo2, nwo3);
                isp = LineHelper.GetIntersectPoint(lineNWOb, linePL, PointType.XZ);

                if (!double.IsNaN(isp.X))
                {
                    CreateNewPointOnPlLine(pl1Line, intersection1);
                    CreateNewPointOnPlLine(pl1Line, isp, cOffsetPhreaticLineBelowSurface, intersection1);
                    CreateNewPointOnPlLine(pl1Line, nwo3Pl);
                    if ((nwo4Pl != null) && (nwo4Pl.X > nwo3Pl.X))
                    {
                        CreateNewPointOnPlLine(pl1Line, nwo4Pl);
                    }
                }
                else
                {
                    string message =
                        "Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.";
                    throw new System.Exception(message);
                }
            }
        }

        /// <summary>
        ///     create new point on geometryPointsString
        /// </summary>
        /// <param name="pl1Line"></param>
        /// <param name="newPoint"></param>
        /// <param name="offsetPhreaticLineBelowSurface"></param>
        /// <param name="zValue"></param>
        private static void CreateNewPointOnPlLine(GeometryPointString pl1Line, GeometryPoint newPoint,
                                                   double offsetPhreaticLineBelowSurface = double.NaN, GeometryPoint zValue = null)
        {
            if (newPoint != null)
            {
                if (double.IsNaN(offsetPhreaticLineBelowSurface))
                {
                    LineHelper.EnsurePointAtX(newPoint.X, pl1Line.Points, GeometryPoint.Precision);
                }
                else
                {
                    LineHelper.EnsurePointAtX(newPoint.X + cOffsetPhreaticLineBelowSurface, pl1Line.Points,
                                              GeometryPoint.Precision);
                }
                if (zValue == null)
                {
                    newPoint.Z = newPoint.Z;
                }
                else
                {
                    newPoint.Z = zValue.Z;
                }
            }
        }

        /// <summary>
        ///     create new line from begin and end point
        /// </summary>
        /// <param name="pointBegin"></param>
        /// <param name="pointEnd"></param>
        /// <returns></returns>
        private static Line CreateNewLine(GeometryPoint pointBegin, GeometryPoint pointEnd)
        {
            return new Line
                {
                    BeginPoint = new GeometryPoint(pointBegin.X, 0, pointBegin.Z),
                    EndPoint = new GeometryPoint(pointEnd.X, 0, pointEnd.Z)
                };
        }

        /// <summary>
        ///     move selected points below surface line
        /// </summary>
        /// <param name="plPointsToBeMoved"></param>
        /// <param name="surfaceline"></param>
        private void MoveSelectedPLLinePointsBelowSurfaceLine(IEnumerable<GeometryPoint> plPointsToBeMoved,
                                                              SurfaceLine2 surfaceline)
        {
            foreach (var plLinePoint in plPointsToBeMoved)
            {
                // Determine which of these points must be moved and move them
                if (surfaceline.Geometry.PositionXzOfPointRelatedToExtrapolatedLine(plLinePoint) !=
                    RelativeXzPosition.BelowGeometricLine)
                {
                    plLinePoint.Z = surfaceline.Geometry.GetZAtX(plLinePoint.X) - cOffsetPhreaticLineBelowSurface;
                }
            }
        }

        /// <summary>
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <param name="lowWaterLevel"></param>
        /// <param name="highWaterLevel"></param>
        private void CreatePhreaticLineSegmentsInsideDikeForLowRiverLevel(GeometryPointString phreaticLine,
                                                                          double lowWaterLevel, Location location)
        {
            GeometryPoint intersectionLowWaterLevelWithDike =
                DetermineIntersectionBetweenWaterLevelAndDike(lowWaterLevel, location.Surfaceline);
            GeometryPoint pointDikeToeAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
            switch (location.PlLineCreationMethod)
            {
                case PlLineCreationMethod.ExpertKnowledgeLinearInDike:
                    const double cPLLineOffsetBelowSurface = 0.1;
                    GeometryPoint intersectionHighWaterLevelWithDike =
                        DetermineIntersectionBetweenWaterLevelAndDike(location.WaterLevelRiver, location.Surfaceline);
                    if (intersectionLowWaterLevelWithDike == null)
                    {
                        // Intersection is supposed to be at toe of dike when no intersection found (i.e. waterlevel below toe of dike)
                        intersectionLowWaterLevelWithDike = new GeometryPoint(pointDikeToeAtRiver.X, 0,
                                                                              pointDikeToeAtRiver.Z);
                    }
                    // At this point both intersectionHighWaterLevelWithDike and intersectionLowWaterLevelWithDike are <> null
                    // or else an exception would have been thrown
                    /*GeometryPoint pointBelowHighLevel = new GeometryPoint(intersectionHighWaterLevelWithDike.X, 0,
                                                                          intersectionHighWaterLevelWithDike.Z -
                                                                          cPLLineOffsetBelowSurface);*/

                    //Add points below surface of dike talud riverside
                    foreach (var point in location.Surfaceline.Geometry.Points.Where(
                        point => point.X > intersectionLowWaterLevelWithDike.X &&
                                 point.X < intersectionHighWaterLevelWithDike.X))
                    {
                        phreaticLine.Points.Add(new GeometryPoint(point.X, 0, point.Z - cPLLineOffsetBelowSurface));
                    }

                    // Add last point (below high riverlevel/dike section
                    //                    phreaticLine.Points.Add(pointBelowHighLevel);

                    break;
                case PlLineCreationMethod.ExpertKnowledgeRrd:
                    //Add points below surface of dike talud riverside until toe of dike riverside
                    foreach (var point in location.Surfaceline.Geometry.Points.Where(
                        point => point.X > intersectionLowWaterLevelWithDike.X &&
                                 point.X <= pointDikeToeAtRiver.X))
                    {
                        if (!IsNonWaterRetainingObjectPoint(point, location.Surfaceline))
                        {
                            phreaticLine.Points.Add(new GeometryPoint(point.X, 0, point.Z - cPLLineOffsetBelowSurface));
                        }
                    }

                    // Add points below crest of dike
                    double offsetDikeTopAtRiver = location.WaterLevelRiver - location.PlLineOffsetBelowDikeTopAtRiver;
                    double offsetDikeTopAtPolder = location.WaterLevelRiver - location.PlLineOffsetBelowDikeTopAtPolder;

                    if (location.Surfaceline.IsDefined(CharacteristicPointType.DikeTopAtRiver))
                    {
                        GeometryPoint pointDikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
                        phreaticLine.Points.Add(new GeometryPoint(pointDikeTopAtRiver.X, 0, offsetDikeTopAtRiver));
                    }

                    if (location.Surfaceline.IsDefined(CharacteristicPointType.DikeTopAtPolder))
                    {
                        GeometryPoint pointDikeTopAtPolder =
                            location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
                        phreaticLine.Points.Add(new GeometryPoint(pointDikeTopAtPolder.X, 0, offsetDikeTopAtPolder));
                    }
                    break;
            }
        }

        /// <summary>
        ///     Create the phreatic line segments inside the dike
        /// </summary>
        /// <param name="phreaticLine"></param>
        private void CreatePhreaticLineSegmentsInsideDikeForHighRiverLevel(GeometryPointString phreaticLine, Location location)
        {
            /*            double offsetDikeTopAtRiver = location.WaterLevelRiverAverage - PLLineOffsetBelowDikeTopAtRiver;
                        double offsetDikeTopAtPolder = location.WaterLevelRiverAverage - PLLineOffsetBelowDikeTopAtPolder;*/
            double offsetDikeTopAtRiver = Math.Round(location.WaterLevelRiver - location.PlLineOffsetBelowDikeTopAtRiver, cRoundDecimals);
            double offsetDikeTopAtPolder = Math.Round(location.WaterLevelRiver - location.PlLineOffsetBelowDikeTopAtPolder, cRoundDecimals);

            GeometryPoint point;

            // points created here are in dike
            switch (location.PlLineCreationMethod)
            {
                case PlLineCreationMethod.ExpertKnowledgeLinearInDike:
                    // No extra points below top of dike at river and top of dike at polder
                    break;
                case PlLineCreationMethod.ExpertKnowledgeRrd:
                    if (location.Surfaceline.IsDefined(CharacteristicPointType.DikeTopAtRiver))
                    {
                        point = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
                        phreaticLine.Points.Add(new GeometryPoint(point.X, 0, offsetDikeTopAtRiver));
                    }

                    if (location.Surfaceline.IsDefined(CharacteristicPointType.DikeTopAtPolder))
                    {
                        point = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
                        phreaticLine.Points.Add(new GeometryPoint(point.X, 0, offsetDikeTopAtPolder));
                    }
                    break;
            }
        }

        /// <summary>
        /// Create the phreatic line segments inside the dike for Phreatic line generated for Ringtoets WTI 2017
        /// </summary>
        /// <param name="phreaticLine">The current phreatic line, in construction</param>
        /// <param name="location">The hydraulic conditions of the waternet creator</param>
        private void CreatePhreaticLineSegmentsInsideDikeForRingtoetsWti2017(GeometryPointString phreaticLine, Location location)
        {
            bool isFinished = false;
            // Shortcuts
            GeometryPoint intersectionPolderLevelWithDike = DetermineIntersectionBetweenPolderLevelAndDike(location.WaterLevelPolder, location.Surfaceline);
            double offsetBelowPointB = 0.0;
            double offsetBelowDikeTopAtPolder = 0.0;
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeTopAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            GeometryPoint dikeToeAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint pointA = phreaticLine.Points.Last();
            double maxXCoordinateSurface = location.Surfaceline.Geometry.Points.Max(x => x.X);

            // Determine the offset values
            if (!location.UseDefaultOffsets)
            {
                offsetBelowPointB = location.PlLineOffsetBelowPointBRingtoetsWti2017;
                offsetBelowDikeTopAtPolder = location.PlLineOffsetBelowDikeTopAtPolder;
            }
            else
            {
                if (location.DikeSoilScenario == DikeSoilScenario.ClayDikeOnClay || location.DikeSoilScenario == DikeSoilScenario.ClayDikeOnSand)
                {
                    offsetBelowPointB = 1.0;
                    offsetBelowDikeTopAtPolder = 1.5;
                }
                if (location.DikeSoilScenario == DikeSoilScenario.SandDikeOnClay)
                {
                    offsetBelowPointB = 0.5 * (pointA.Z - dikeToeAtPolder.Z);
                }
            }

            // Determine X coordinate of point B
            double xPointB;
            if (location.DikeSoilScenario == DikeSoilScenario.SandDikeOnClay)
            {
                xPointB = Math.Min(phreaticLine.Points.Last().X + 0.001, dikeTopAtPolder.X);
            }
            else
            {
                xPointB = Math.Min(phreaticLine.Points.Last().X + 1.0, dikeTopAtPolder.X);
            }

            // Determine the minimum value for point B (linear interpolation between both minimum values at Dike Top)
            double minValueForPointB = (location.MinimumLevelPhreaticLineAtDikeTopPolder - location.MinimumLevelPhreaticLineAtDikeTopRiver)
                / (dikeTopAtPolder.X - dikeTopAtRiver.X) * (xPointB - dikeTopAtRiver.X) + location.MinimumLevelPhreaticLineAtDikeTopRiver;

            // Determine the phreatic levels at points B and DikeTopPolder
            double zPointB = Math.Round(Math.Min(location.Surfaceline.Geometry.GetZAtX(xPointB) - cOffsetPhreaticLineBelowSurface, 
                Math.Max(location.WaterLevelRiver - offsetBelowPointB, minValueForPointB)), cRoundDecimals);
            double waterLevelBelowDikeTopAtPolder = Math.Max(Math.Round(location.WaterLevelRiver - offsetBelowDikeTopAtPolder, cRoundDecimals),
                                                             location.MinimumLevelPhreaticLineAtDikeTopPolder);
            
            // Both created points are now added to the phreatic line. However:
            // - For Sand dike on sand (with default offset), no points are added below both dike top yet;
            // - For Sand dike on clay (with default offset), only 1 point below dike top at river side is added.
            // This is done afterwards, in method CreatePLLine1ByExpertKnowledge, because the position of point E has to be known.
            if (location.Surfaceline.IsDefined(CharacteristicPointType.DikeTopAtRiver))
            {
                // If the water level at polder if above the phreatic line point at the top of the dike at river side (point B), then the rest of the phreatic line
                // is an horizontal line at water polder, and starting at intersection polder level/dike
                if (zPointB < location.WaterLevelPolder)
                {
                    if (intersectionPolderLevelWithDike != null)
                    {
                        phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                        phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                        isFinished = true;
                    }
                    else
                    {
                        if (!(location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand && location.UseDefaultOffsets))
                        {
                            phreaticLine.Points.Add(new GeometryPoint(xPointB, 0, location.WaterLevelPolder));
                        }
                    }
                }
                else
                {
                    if (!(location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand && location.UseDefaultOffsets))
                    {
                        phreaticLine.Points.Add(new GeometryPoint(xPointB, 0, zPointB));
                    }
                }
            }

            if (location.Surfaceline.IsDefined(CharacteristicPointType.DikeTopAtPolder) && !isFinished)
            {
                // If the water level at polder if above the phreatic line point at the top of the dike at river side (point B), then the rest of the phreatic line
                // is an horizontal line at water polder, and starting at intersection polder level/dike
                if (waterLevelBelowDikeTopAtPolder < location.WaterLevelPolder)
                {
                    if (intersectionPolderLevelWithDike != null)
                    {
                        phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                        phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                    }
                    else
                    {
                        phreaticLine.Points.Add(new GeometryPoint(dikeTopAtPolder.X, 0, location.WaterLevelPolder));
                    }
                }
                else
                {
                    if (!(location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand && location.UseDefaultOffsets))
                    {
                        if (!(location.DikeSoilScenario == DikeSoilScenario.SandDikeOnClay && location.UseDefaultOffsets))
                        {
                            phreaticLine.Points.Add(new GeometryPoint(dikeTopAtPolder.X, 0, waterLevelBelowDikeTopAtPolder));
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Create the phreatic line segments inside the dike for Phreatic line generated for Ringtoets WTI 2017
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <param name="location"> </param>
        private void CreatePhreaticLineSegmentsInShoulderAndPolderForRingtoetsWti2017(GeometryPointString phreaticLine, Location location)
        {
            bool isFinished = false;
            // Shortcuts
            GeometryPoint intersectionPolderLevelWithDike = DetermineIntersectionBetweenPolderLevelAndDike(location.WaterLevelPolder, location.Surfaceline);
            GeometryPoint shoulderBaseInsidePoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside);
            GeometryPoint dikeToeAtPolderPoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint ditchDikeSidePoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
            bool isDitchPresent = (ditchDikeSidePoint != null && !double.IsNaN(ditchDikeSidePoint.X));
            double maxXCoordinateSurface = location.Surfaceline.Geometry.Points.Max(x => x.X);
            GeometryPoint intersectionDitchWithWaterPolder = DetermineIntersectionBetweenPolderLevelAndDitch(location, isDitchPresent);

            // If the water level at polder if above the phreatic line point at the top of the dike (point C), then stop
            if (Math.Abs(phreaticLine.Points.Last().X - maxXCoordinateSurface) < 0.00001)
            {
                isFinished = true;
            }

            // If a drainage construction is present, the phreatic line goes directly to the middle of it (shoulder offset not used)
            if (!isFinished)
            {
                if ((location.DikeSoilScenario == DikeSoilScenario.SandDikeOnClay ||
                     location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand) && location.DrainageConstructionPresent)
                {
                    ThrowWhenDrainageConstructionNotInsideDike(location, location.Surfaceline);
                    phreaticLine.Points.Add(new GeometryPoint(location.XCoordMiddleDrainageConstruction, 0,
                                                              location.ZCoordMiddleDrainageConstruction));
                    // If a ditch is present, go directly to intersection point with water polder level
                    if (isDitchPresent && !double.IsNaN(intersectionDitchWithWaterPolder.X))
                    {
                        phreaticLine.Points.Add(intersectionDitchWithWaterPolder);
                        phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                        isFinished = true;
                    }
                }
                else
                {
                    // A point is created below shoulder point only for user-defined offset, and only if a shoulder is present
                    if (location.Surfaceline.IsDefined(CharacteristicPointType.ShoulderBaseInside) && !location.UseDefaultOffsets)
                    {
                        // determine offset value below the shoulder: a negative offset is not allowed, in this case the method uses 1 cm
                        double offsetAtShouder = Math.Max(cOffsetPhreaticLineBelowSurface, location.PlLineOffsetBelowShoulderBaseInside);
                        double zLevelAtShoulder = shoulderBaseInsidePoint.Z - offsetAtShouder;
                        // automatic correction of the phreatic line to avoid that it goes up 
                        zLevelAtShoulder = Math.Min(phreaticLine.Points.Last().Z, zLevelAtShoulder);
                        // If the water level at polder if above the phreatic line point at shoulder, then the rest of the phreatic line
                        // is an horizontal line at water polder, and starting at intersection polder level/dike
                        if (zLevelAtShoulder < location.WaterLevelPolder)
                        {
                            phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                            phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                            isFinished = true;
                        }
                        else
                        {
                            phreaticLine.Points.Add(new GeometryPoint(shoulderBaseInsidePoint.X, 0, zLevelAtShoulder));
                        }
                    }
                    // For sand dikes, an extra point (D1) is created 0.25*h above dike toe, where h is the distance between
                    // point A (water level at river side) and point E (dike toe at polder side)
                    bool negativeUserDefinedOffset = (!location.UseDefaultOffsets && (location.PlLineOffsetBelowDikeToeAtPolder < 0));
                    bool sandDikesWithDefaultOffsets = (location.DikeSoilScenario == DikeSoilScenario.SandDikeOnClay || location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand)
                        && location.UseDefaultOffsets;
                    if ((sandDikesWithDefaultOffsets || negativeUserDefinedOffset) && !isFinished)
                    {
                        double levelPointA = phreaticLine.Points.First().Z;
                        double levelPointD1;
                        if (location.UseDefaultOffsets)
                        {
                            levelPointD1 = dikeToeAtPolderPoint.Z + 0.25 * (levelPointA - dikeToeAtPolderPoint.Z);
                        }
                        else
                        {
                            levelPointD1 = dikeToeAtPolderPoint.Z - location.PlLineOffsetBelowDikeToeAtPolder;
                        }
                        GeometryPoint pointD1 = DetermineIntersectionBetweenPolderLevelAndDike(levelPointD1, location.Surfaceline);
                        if (levelPointD1 < location.WaterLevelPolder)
                        {
                            phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                            phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                            isFinished = true;
                        }
                        else
                        {
                            if (pointD1 != null)
                            {
                                phreaticLine.Points.Add(new GeometryPoint(pointD1.X, 0, pointD1.Z - cOffsetPhreaticLineBelowSurface));
                            }
                        }
                    }
                }
            }

            // If the water level at polder if above the phreatic line point at shoulder, then the rest of the phreatic line
            // is an horizontal line at water polder, and starting at intersection polder level/dike
            if (dikeToeAtPolderPoint.Z < location.WaterLevelPolder && !isFinished)
            {
                phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                isFinished = true;
            }

            // Add point E (dike toe at polder side)
            if (!Double.IsNaN(dikeToeAtPolderPoint.X) && !isFinished)
            {
                // determine offset value below the dike toe polder: if a negative offset is used (i.e. point D1 is created),
                // point E is situated 1 cm below dike toe
                double offsetBelowDikeToe = 0.0;
                if (!location.UseDefaultOffsets)
                {
                    offsetBelowDikeToe = Math.Max(cOffsetPhreaticLineBelowSurface, location.PlLineOffsetBelowDikeToeAtPolder);
                }

                double zLevelAtDikeToe = dikeToeAtPolderPoint.Z - offsetBelowDikeToe;
                // automatic correction of the phreatic line to avoid that it goes up between point E and F, if point F (intersection with ditch) exists
                if (isDitchPresent && !double.IsNaN(intersectionDitchWithWaterPolder.X))
                {
                    zLevelAtDikeToe = Math.Max(intersectionDitchWithWaterPolder.Z, zLevelAtDikeToe);
                }
                else
                {
                    // if no ditch present, or no intersection with ditch
                    offsetBelowDikeToe = Math.Max(cOffsetPhreaticLineBelowSurface, offsetBelowDikeToe);
                    zLevelAtDikeToe = Math.Max(location.WaterLevelPolder, dikeToeAtPolderPoint.Z - offsetBelowDikeToe);
                }
                phreaticLine.Points.Add(new GeometryPoint(dikeToeAtPolderPoint.X, 0, zLevelAtDikeToe));
            }

            // If point F (intersection point between the ditch slope at dike side and the level of the water polder) exists, add it to phreatic line
            if (isDitchPresent && !double.IsNaN(intersectionDitchWithWaterPolder.X) && !isFinished)
            {
                //Todo #Bka see MAC-242: This point should not just be added. First you should see if there are other points 
                //on the surfaceline between the last added phreatic point and this intersectionDitchWithWaterPolder. If there
                // are such points, see if the intended line does not fall above these points. If so, then add point(s) to the
                // phreatic line at that X with the give offset below the point(s).
                // It is in fact easier to just add a phreatic point for all points that fall between the last added phreatic 
                // point and this intersectionDitchWithWaterPolder (with offset of course) but I do not know if that is intended.
                // For that check with Alexander van Duinen.
                phreaticLine.Points.Add(intersectionDitchWithWaterPolder);
            }
            // If no ditch is present and if the phreatic level at dike toe (i.e. point E = dike toe - offset)
            // is above the water level at polder, then the Phreatic Line follows the surface line minus offset (= 1 cm).
            else
            {
                if (!isFinished && (phreaticLine.Points[phreaticLine.Points.Count - 1].Z > location.WaterLevelPolder))
                {
                    AddPhreaticLineAlongSurfaceLevel(phreaticLine, location);
                }
            }

            // Validate if end point surface is reached
            if (Math.Abs(phreaticLine.Points.Last().X - maxXCoordinateSurface) > 0.0001 && !isFinished)
            {
                var endPoint = new GeometryPoint(maxXCoordinateSurface, 0, phreaticLine.Points.Last().Z);
                phreaticLine.Points.Add(endPoint);
            }
        }

        /// <summary>
        /// If a ditch is present, determine point F = the intersection point between the ditch slope at dike side and the level of the water polder.
        /// If no point is found, return NaN.
        /// </summary>
        /// <param name="location">Class containing hydraulic data's over the waternet creator</param>
        /// <param name="isDitchPresent">is a ditch present ?</param>
        /// <returns> The intersection point between the ditch slope at dike side and the level of the water polder. If no point is found, return NaN.</returns>
        private static GeometryPoint DetermineIntersectionBetweenPolderLevelAndDitch(Location location, bool isDitchPresent)
        {
            GeometryPoint ditchDikeSidePoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
            GeometryPoint bottomDitchPolderSidePoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchPolderSide);
            double minXCoordinateSurface = location.Surfaceline.Geometry.Points.Min(x => x.X);
            double maxXCoordinateSurface = location.Surfaceline.Geometry.Points.Max(x => x.X);
            var intersectionDitchWithWaterPolder = new GeometryPoint(double.NaN, double.NaN, double.NaN);
            if (isDitchPresent)
            {
                var waterLevelPolderLine = new Line();
                var firstPoint = new GeometryPoint(minXCoordinateSurface, 0, location.WaterLevelPolder);
                var lastPoint = new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder);
                waterLevelPolderLine.SetBeginAndEndPoints(firstPoint, lastPoint);

                GeometryPointString ditchPart = location.Surfaceline.Geometry.GetPart(ditchDikeSidePoint.X, bottomDitchPolderSidePoint.X);
                for (int i = 0; i < ditchPart.Count - 1; i++)
                {
                    var lineDitchDikeSide = new Line(ditchPart[i], ditchPart[i + 1]);
                    intersectionDitchWithWaterPolder = LineHelper.GetIntersectPoint(lineDitchDikeSide, waterLevelPolderLine, PointType.XZ);
                    if (!double.IsNaN(intersectionDitchWithWaterPolder.X))
                    {
                        break;
                    }
                }
            }
            return intersectionDitchWithWaterPolder;
        }

        /// <summary>
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <param name="location"> </param>
        private void CreatePhreaticLineSegmentsInShoulderAndPolder(GeometryPointString phreaticLine, Location location)
        {
            bool isFinished = false;
            GeometryPoint intersectionPolderLevelWithDike =
                DetermineIntersectionBetweenPolderLevelAndDike(location.WaterLevelPolder, location.Surfaceline);

            double maxXCoordinateSurface = location.Surfaceline.Geometry.Points.Max(x => x.X);
            if (location.PlLineCreationMethod != PlLineCreationMethod.ExpertKnowledgeLinearInDike)
            {
                // Points created below are points starting at shoulder point to the right
                if (location.Surfaceline.IsDefined(CharacteristicPointType.ShoulderBaseInside))
                {
                    GeometryPoint shoulderBaseInsidePoint =
                        location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside);
                    double zLevel = Math.Min(phreaticLine.Points.Last().Z,
                                             shoulderBaseInsidePoint.Z -
                                             Math.Max(cOffsetPhreaticLineBelowSurface, location.PlLineOffsetBelowShoulderBaseInside));
                    if (zLevel < location.WaterLevelPolder)
                    {
                        // if polderlevel above point rest of freatic line is polderlevel starting with intersection polderlevel dike
                        phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                        phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                        isFinished = true;
                    }
                    else
                    {
                        phreaticLine.Points.Add(new GeometryPoint(shoulderBaseInsidePoint.X, 0, zLevel));
                    }
                }
            }

            GeometryPoint dikeToeAtPolderPoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
            GeometryPoint ditchDikeSidePoint = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
            if (dikeToeAtPolderPoint != null && !Double.IsNaN(dikeToeAtPolderPoint.X) && (!isFinished))
            {
                double zLevel = Math.Min(phreaticLine.Points.Last().Z,
                                         Math.Round(dikeToeAtPolderPoint.Z - Math.Max(cOffsetPhreaticLineBelowSurface, location.PlLineOffsetBelowDikeToeAtPolder), cRoundDecimals));
                if (zLevel < location.WaterLevelPolder)
                {
                    // if polderlevel above point rest of freatic line is polderlevel starting with intersection polderlevel dike
                    if (!(intersectionPolderLevelWithDike == null))
                    {
                        phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                    }
                    else
                    {
                        phreaticLine.Points.Add(new GeometryPoint(dikeToeAtPolderPoint.X, 0, location.WaterLevelPolder));
                    }
                    phreaticLine.Points.Add(new GeometryPoint(maxXCoordinateSurface, 0, location.WaterLevelPolder));
                    isFinished = true;
                }
                else
                {
                    if (ditchDikeSidePoint != null && !Double.IsNaN(ditchDikeSidePoint.X))
                    {
                        if (ditchDikeSidePoint.LocationEquals(dikeToeAtPolderPoint))
                        {
                            // If DikeToeAtPolder is same as DitchDikeSide pl1 should always go to polderlevel at this point
                            zLevel = location.WaterLevelPolder;
                        }
                    }

                    zLevel = Math.Max(zLevel, location.WaterLevelPolder);
                    // Add point if it lies left of intersection of polderlevel with dike
                    if ((intersectionPolderLevelWithDike == null) ||
                        (intersectionPolderLevelWithDike.X > dikeToeAtPolderPoint.X))
                    {
                        phreaticLine.Points.Add(new GeometryPoint(dikeToeAtPolderPoint.X, 0, zLevel));
                    }
                }
            }

            if (!isFinished)
            {
                if (intersectionPolderLevelWithDike != null)
                {
                    phreaticLine.Points.Add(intersectionPolderLevelWithDike);
                }

                bool isDitchPresent = (ditchDikeSidePoint != null && !double.IsNaN(ditchDikeSidePoint.X));
                bool isNonWaterRetainingOjectPresent = ((location.NWOPhreaticAdaption != PhreaticAdaptionType.None) &&
                                                        (location.Surfaceline.IsDefined(
                                                            CharacteristicPointType.NonWaterRetainingObjectPoint1)));
                bool adjustDitch = false;
                // Handle making the waterlevel horizontal in the NWO at the Polderside when needed (For Riverside see AdaptPL1ForNonWaterRetainingObject).
                if (isNonWaterRetainingOjectPresent && (location.NWOPhreaticAdaption != PhreaticAdaptionType.Fill) &&
                    (location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X >=
                     location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X))
                {
                    // if there is a ditch and it is to the left of the NWO, then only the ditch needs to be adjusted and the NWO will be correct automatically.
                    // if there is a ditch but it is to the right of the NWO, the NWO needs adjusting and the Ditch will be correct automatically.
                    // if there is no ditch then the NWO needs adjusting.
                    if (isDitchPresent)
                    {
                        adjustDitch = (location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X >=
                                       location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X);
                    }
                    if (!adjustDitch)
                    {
                        GeometryPoint nw1 = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1);
                        AdjustForDitchAndOrNonWaterRetainingObjectatPolderSide(phreaticLine, nw1.X, location);
                    }
                }
                else
                {
                    // No (relevant) NWO so enable handling of ditch.
                    // If there is a ditch but there is also a NWO to the right of it at polder side, then do not adjust the ditch. Do in all other cases.
                    // First see if there is a NWO.
                    adjustDitch = !isNonWaterRetainingOjectPresent;
                    if (!adjustDitch)
                    {
                        // there is a NWO, check the position
                        adjustDitch =
                            !((isDitchPresent) &&
                              (location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X >=
                               location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X) &&
                              (location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).X <=
                               location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X));
                    }
                }
                // if the NWO is not there or irrelevant and there is a ditch, then adjust it.
                if ((adjustDitch) && (isDitchPresent))
                {
                    AdjustForDitchAndOrNonWaterRetainingObjectatPolderSide(phreaticLine, ditchDikeSidePoint.X, location);
                }
                else
                {
                    // if no ditch then the PL1 will continue horizontally until the end of the profile
                    // another choice could be to let the PL1 go to polderlevel at the end of the profile, but that could be too optimistic for uplift.
                    // Discussion about this was done between Vastenburg, Knoeff and The.
                    // After a renewed discussion with Vastenburg, Van der Zwan and Bka, it has been decided that the PL1 should follow the 
                    // surfaceline (with offset) until either end of surface line or polder level. Note: this is only needed when the phreatic level
                    // at dike toe is above polder level.
                    if ((!isNonWaterRetainingOjectPresent) && (!isDitchPresent))
                    {
                        if (phreaticLine.Points[phreaticLine.Points.Count - 1].Z > location.WaterLevelPolder)
                        {
                            AddPhreaticLineAlongSurfaceLevel(phreaticLine, location);
                        }
                    }
                }

                //Validate if endpoint surface has reached
                if (phreaticLine.Points.Last().X != maxXCoordinateSurface)
                {
                    var endPoint = new GeometryPoint(maxXCoordinateSurface, 0, phreaticLine.Points.Last().Z);
                    phreaticLine.Points.Add(endPoint);
                }
            }
        }

        /// <summary>
        /// Add phreatic point at offset (i.e. point E = DikeToeAtPolder - offset) below every surface line point as long as depth > polder level. Else determine the
        /// proper position of the point at polder level (intersection) and stop.
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <param name="location"></param>
        private void AddPhreaticLineAlongSurfaceLevel(GeometryPointString phreaticLine, Location location)
        {
            GeometryPointString relevantPart =
                location.Surfaceline.Geometry.GetPart(
                    location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, location.Surfaceline.Geometry.GetMaxX());

            bool intersected = false;
            double offset;
            if (location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017 && location.UseDefaultOffsets)
            {
                offset = 0;
            }
            else
            {
                offset = location.PlLineOffsetBelowDikeToeAtPolder;
            }
            offset = Math.Max(cOffsetPhreaticLineBelowSurface, offset);

            for (int i = 0; i < relevantPart.Count; i++)
            {
                double z = Math.Max(location.WaterLevelPolder, relevantPart.Points[i].Z - offset);
                double x = relevantPart.Points[i].X;

                if (location.WaterLevelPolder > relevantPart[i].Z - offset)
                {
                    // Determine intersection between would be phreatic segment and polderlevel. Add that as next point.
                    var waterLevelPolderLine = new Line
                        {
                            BeginPoint = new GeometryPoint(location.Surfaceline.Geometry.GetMinX(), 0, location.WaterLevelPolder),
                            EndPoint = new GeometryPoint(location.Surfaceline.Geometry.GetMaxX(), 0, location.WaterLevelPolder)
                        };

                    var slopeLine = new Line
                        {
                            BeginPoint = phreaticLine.Points[phreaticLine.Points.Count - 1],
                            EndPoint = new GeometryPoint(relevantPart.Points[i].X, 0, relevantPart.Points[i].Z - offset)
                        };

                    GeometryPoint intersectionPoint = LineHelper.GetIntersectPoint(waterLevelPolderLine, slopeLine, PointType.XZ);
                    if (!double.IsNaN(intersectionPoint.X))
                    {
                        x = intersectionPoint.X;
                    }
                    intersected = true;
                }
                var point = new GeometryPoint(x, 0, z);
                phreaticLine.Points.Add(point);
                if (intersected)
                {
                    break;
                }
            }
        }

        private void AdjustForDitchAndOrNonWaterRetainingObjectatPolderSide(GeometryPointString phreaticLine, double x, Location location)
        {
            const double maxDouble = 99999.999;
            var partialLine = new Line();
            //#bka: hier niet langer ook starten met waterlevel als waterlevel onder bottomditch zit!
            partialLine.SetBeginAndEndPoints(new GeometryPoint(phreaticLine.Points[0].X, 0, location.WaterLevelPolder),
                                             new GeometryPoint(maxDouble, 0, location.WaterLevelPolder));
            AddIntersectionDitchDikeSegmentPolderLevelToPhreatic(phreaticLine, x, partialLine, location.Surfaceline);
            AddIntersectionDitchPolderSegmentPolderLevelToPhreatic(phreaticLine, partialLine, location.Surfaceline);
        }

        /// <summary>
        ///     Intersection between two line segments:
        ///     -Ditchpolder Surfaceline segment
        ///     -Polder waterlevel
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <param name="phreaticPolderPartialLine"></param>
        /// <param name="surfaceline"> </param>
        private void AddIntersectionDitchPolderSegmentPolderLevelToPhreatic(GeometryPointString phreaticLine,
                                                                            Line phreaticPolderPartialLine, SurfaceLine2 surfaceLine)
        {
            if (surfaceLine.IsDefined(CharacteristicPointType.DitchPolderSide))
            {
                GeometryPoint pointDitchPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide);

                GeometryPointString relevantPart = surfaceLine.Geometry.GetPart(surfaceLine.Geometry.GetMinX(), pointDitchPolder.X);

                var lineDitchPolderSide = new Line();
                if (relevantPart.Points.Count > 1)
                {
                    lineDitchPolderSide.BeginPoint = relevantPart.Points[relevantPart.Points.Count - 2];
                    lineDitchPolderSide.EndPoint = relevantPart.Points[relevantPart.Points.Count - 1];

                    GeometryPoint intersectDitchPolderPhreatic = LineHelper.GetIntersectPoint(lineDitchPolderSide, phreaticPolderPartialLine, PointType.XZ);
                    if (!double.IsNaN(intersectDitchPolderPhreatic.X))
                    {
                        phreaticLine.Points.Add(new GeometryPoint(intersectDitchPolderPhreatic.X, 0, intersectDitchPolderPhreatic.Z));
                    }
                    else
                    {
                        phreaticLine.Points.Add(new GeometryPoint(pointDitchPolder.X, 0, phreaticPolderPartialLine.EndPoint.Z));
                    }
                }
                else
                {
                    string message = "DitchPolderSide should be part of a line segment";
                    throw new System.Exception(message);
                }
            }
        }

        /// <summary>
        ///     Intersection between two line segments:
        ///     -DitchDike Surfaceline segment
        ///     -Polder waterlevel
        /// </summary>
        /// <param name="phreatic"></param>
        /// <param name="x">Ditch x coordinate</param>
        /// <param name="partialLine"></param>
        /// <param name="surfaceLine"> </param>
        private void AddIntersectionDitchDikeSegmentPolderLevelToPhreatic(GeometryPointString phreatic, double x, Line partialLine, SurfaceLine2 surfaceLine)
        {
            if (x < surfaceLine.Geometry.GetMaxX())
            {
                GeometryPointString relevantPart = surfaceLine.Geometry.GetPart(x, surfaceLine.Geometry.GetMaxX());

                var lineDitchDikeSide = new Line(relevantPart[0], relevantPart[1]);

                GeometryPoint isp = LineHelper.GetIntersectPoint(lineDitchDikeSide, partialLine, PointType.XZ);
                if (!double.IsNaN(isp.X))
                {
                    phreatic.Points.Add(new GeometryPoint(isp.X, 0, isp.Z));
                }
                else
                {
                    phreatic.Points.Add(new GeometryPoint(x, 0, partialLine.EndPoint.Z));
                }
            }
            else
            {
                string message = "Could not create DikeSementPolderLevel";
                throw new System.Exception(message);
            }
        }

        /// <summary>
        /// Determine if polder level intersects the dike at polder side (between dike top at river and dike toe at polder)
        /// </summary>
        /// <param name="polderLevel"></param>
        /// <param name="surfaceline"></param>
        /// <returns></returns>
        private GeometryPoint DetermineIntersectionBetweenPolderLevelAndDike(double polderLevel, SurfaceLine2 surfaceline)
        {
            var polderlevelLine = new Line();
            var localSurfaceline = new GeometryPointString();
            localSurfaceline.Points.AddRange(surfaceline.Geometry.Points.Where(p => !double.IsNaN(p.X)));

            double startXCoordinate = localSurfaceline.Points.OrderBy(p => p.X).First().X;
            double endxCoordinate = localSurfaceline.Points.OrderBy(p => p.X).Last().X;
            polderlevelLine.SetBeginAndEndPoints(new GeometryPoint(startXCoordinate, 0, polderLevel),
                                                 new GeometryPoint(endxCoordinate, 0, polderLevel));

            ThrowWhenWaterLevelPolderAboveDikeTopAtPolder(polderLevel, surfaceline);

            GeometryPointString relevantPart =
                surfaceline.Geometry.GetPart(
                    surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X,
                    surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X);

            for (int index = 0; index < relevantPart.Points.Count - 1; index++)
            {
                var surfaceLineSegment = new Line(relevantPart[index], relevantPart[index + 1]);

                GeometryPoint intersectPoint = LineHelper.GetIntersectPoint(surfaceLineSegment, polderlevelLine, PointType.XZ);
                if (!double.IsNaN(intersectPoint.X))
                {
                    return new GeometryPoint(intersectPoint.X, 0, intersectPoint.Z);
                }
            }
            return null;
        }

        /// <summary>
        /// Detects whether the phreatic line is above the surface line between entry point in dike and polder.
        /// If the phreatic line is above the surface, then the phreatic line should follow the surface
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <returns></returns>
        private void ValidatePhreaticLineBelowDike(GeometryPointString phreaticLine, Location location)
        {
            GeometryPointString relevantPart =
                location.Surfaceline.Geometry.GetPart(phreaticLine[1].X + 0.0001, location.Surfaceline.Geometry.GetMaxX());

            bool stopIteration = false;
            // This code can go into an endless loop
            // Enter failsave to raise exception when this occurs
            int cMaxIterations = Math.Max(100, phreaticLine.Points.Count * location.Surfaceline.Geometry.Count);
            int iterationCount = 0;
            while (stopIteration == false)
            {
                iterationCount++;
                bool foundIntersect = false;
                for (int index = 0; index < phreaticLine.Points.Count - 1; index++)
                {
                    var phreaticLineSegment = new Line(phreaticLine.Points[index], phreaticLine.Points[index + 1]);

                    for (int surfacePointIndex = 0; surfacePointIndex < relevantPart.Count - 1; surfacePointIndex++)
                    {
                        var surfaceLineSegment = new Line(relevantPart.Points[surfacePointIndex], relevantPart.Points[surfacePointIndex + 1]);

                        GeometryPoint intersectPoint = LineHelper.GetIntersectPoint(phreaticLineSegment, surfaceLineSegment, PointType.XZ);

                        if (!double.IsNaN(intersectPoint.X))
                        {
                            // Only add intersectPoint if it is higher than waterlevelPolder, because the intersection could be caused 
                            // by the fact that the waterlevel is higher than DikeToeAtPolder
                            if (intersectPoint.Z > location.WaterLevelPolder + cOffsetPhreaticLineBelowSurface)
                            {
                                var newPLLinePoint = new GeometryPoint(intersectPoint.X, 0,
                                                                       intersectPoint.Z -
                                                                       cOffsetPhreaticLineBelowSurface);
                                // Only add new point if it is more to the right than the first point of the phreatic line segment when
                                if (newPLLinePoint.X > phreaticLineSegment.BeginPoint.X)
                                {
                                    phreaticLine.Points.Insert(index + 1, newPLLinePoint);
                                    if (Equal(relevantPart.Points[surfacePointIndex + 1].X, intersectPoint.X))
                                    {
                                        // phreatic point and surface line point are positioned above each other, so replace the phreatic point with the intersect point
                                        phreaticLine.Points[index + 2] =
                                            new GeometryPoint(relevantPart[surfacePointIndex + 1].X, 0,
                                                              relevantPart.Points[surfacePointIndex + 1].Z -
                                                              cOffsetPhreaticLineBelowSurface);
                                    }
                                    if (Equal(relevantPart.Points[surfacePointIndex].X, intersectPoint.X))
                                    {
                                        // phreatic point and surface line point are positioned above each other, so replace the phreatic point with the intersect point
                                        phreaticLine.Points[index + 2] =
                                            new GeometryPoint(relevantPart[surfacePointIndex].X, 0,
                                                              relevantPart.Points[surfacePointIndex].Z -
                                                              cOffsetPhreaticLineBelowSurface);
                                    }
                                    else
                                    {
                                        var newNextPLLinePoint =
                                            new GeometryPoint(relevantPart.Points[surfacePointIndex + 1].X, 0,
                                                              relevantPart.Points[surfacePointIndex + 1].Z -
                                                              cOffsetPhreaticLineBelowSurface);
                                        if (newNextPLLinePoint.X <= phreaticLine.Points[index + 2].X)
                                        {
                                            if (Equal(newNextPLLinePoint.X,
                                                      phreaticLine.Points[index + 2].X))
                                            {
                                                // If phreatic point already exist on this x-coordinate just replace it with the new point
                                                phreaticLine.Points[index + 2] = newNextPLLinePoint;
                                            }
                                            else
                                            {
                                                phreaticLine.Points.Insert(index + 2, newNextPLLinePoint);
                                            }
                                        }
                                    }
                                    foundIntersect = true;
                                }
                                break;
                            }
                        }
                    }
                }
                if (foundIntersect == false)
                {
                    stopIteration = true;
                }
                if (iterationCount > cMaxIterations)
                {
                    string message = "PLLinesCreator.ValidatePhreaticLineBelowDike() cannot succeed";
                    throw new System.Exception(message);
                }
            }
        }
        

        /// <summary>
        /// Check if the phreatic line is above the water level at polder. If it is below, correct it.
        /// </summary>
        /// <param name="phreaticLine"></param>
        /// <param name="location"></param>
        private void ValidatePhreaticAboveWaterLevelPolder(GeometryPointString phreaticLine, Location location)
        {
            foreach (var phreaticPoint in phreaticLine.Points)
            {
                if (phreaticPoint.Z < location.WaterLevelPolder)
                {
                    phreaticPoint.Z = location.WaterLevelPolder;
                }
            }
        }

        /// <summary>
        /// Check if water level is below dike
        /// </summary>
        /// <param name="waterLevelRiver"></param>
        /// <param name="surfaceLine"></param>
        private void ThrowWhenWaterLevelAboveDike(double waterLevelRiver, SurfaceLine2 surfaceLine)
        {
            if (waterLevelRiver > surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).Z)
            {
                string message =
                    String.Format(
                        "Phreatic water level should have an intersection with the dike at least once (phreatic line higher ({0:F2} m) than surface line ({1:F2}))",
                        waterLevelRiver, surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).Z);
                throw new System.Exception(message);
            }
        }

        /// <summary>
        /// Check if water level at polder is above dike top at polder
        /// </summary>
        private void ThrowWhenWaterLevelPolderAboveDikeTopAtPolder(double polderLevel, SurfaceLine2 surfaceLine)
        {
            if (polderLevel > surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).Z)
            {
                string message =
                    String.Format(
                        "Waterlevel ({0:0.00}) in polder higher than dike top at polder ({1:0.00}) in surface line {2}",
                        polderLevel, GetDikeToeInward(surfaceLine).Z, surfaceLine.Name);
                throw new System.Exception(message);
            }
        }

        /// <summary>
        /// Check if the middle of the drainage construction is situated between DikeTopAtPolder and DikeToeAtPolder.
        /// Check also that the middle of the drainage construction is inside the dike, not in the air.
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        private void ThrowWhenDrainageConstructionNotInsideDike(Location location, SurfaceLine2 surfaceLine)
        {
            double xCoordDrain = location.XCoordMiddleDrainageConstruction;
            double zCoordDrain = location.ZCoordMiddleDrainageConstruction;
            double zTopDike = surfaceLine.Geometry.GetZAtX(xCoordDrain);
            double xDikeToeAtPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
            double xDikeTopAtPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).X;
            if ((xCoordDrain > xDikeToeAtPolder) || (xCoordDrain < xDikeTopAtPolder))
            {
                string message =
                    String.Format(
                        "The middle of the drainage construction ({0:0.00}) must be situated between the dike top ({1:0.00}) and the dike toe at polder side ({2:0.00}) " +
                        "in surface line {3}", xCoordDrain, xDikeTopAtPolder, xDikeToeAtPolder, surfaceLine.Name);
                throw new System.Exception(message);
            }
            if (zCoordDrain > zTopDike)
            {
                string message =
                    String.Format(
                        "The drainage construction ({0:0.00}) must be situated below surface line {1}", zCoordDrain, surfaceLine.Name);
                throw new System.Exception(message);
            }
        }
    }
}