// Copyright (C) Stichting Deltares 2016. All rights reserved.
//
// This file is part of the Macro Stability kernel.
//
// The Macro Stability kernel is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// All names, logos, and references to "Deltares" are registered trademarks of
// Stichting Deltares and remain full property of Stichting Deltares at all times.
// All rights reserved.

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

namespace Deltares.Geotechnics.WaternetCreator
{
    /// <summary>
    /// This class creates 4 PL-lines (PL1, PL2, PL3 and PL4) for method Ringtoets WTI 2017.
    /// Refer to the memo of Alexander van Duinen (1209434-012-GEO-0002-m-Schematisering waterspanningen in WTI 2017 (Ringtoets).docx) for more information:
    /// https://repos.deltares.nl/repos/FailureMechanisms/FailureMechanisms/DikesMacroStability/trunk/doc/Work
    /// </summary>
    public class PlLinesCreatorWTI
    {
        private const double cUpliftFactorEquilibrium = 1.0;
        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 phreaticLine;
        private GeometryPointString Pl2;

        /// <summary>
        /// Initializes a new instance of the <see cref="PlLinesCreator"/> class.
        /// </summary>
        public PlLinesCreatorWTI()
        {
            InitPLLinesCreator();
        }

        /// <summary>
        /// Gets or sets a value indicating whether oven dry unit weight is used.
        /// </summary>
        /// <value>
        /// <c>true</c> if oven dry unit weight is used; otherwise, <c>false</c>.
        /// </value>
        public bool IsUseOvenDryUnitWeight { get; set; }

        /// <summary>
        /// Gets or sets the type of the soil geometry.
        /// </summary>
        /// <value>
        /// The type of the soil geometry.
        /// </value>
        public SoilGeometryType SoilGeometryType { get; set; }

        /// <summary>
        /// Gets or sets the dike embankment material.
        /// </summary>
        /// <value>
        /// The dike embankment material.
        /// </value>
        public Soil DikeEmbankmentMaterial { get; set; }

        /// <summary>
        /// Gets or sets the head in PL-line 2.
        /// </summary>
        /// <value>
        /// The head in PL-line 2.
        /// </value>
        public double? HeadInPLLine2 { get; set; }

        /// <summary>
        /// Gets or sets the unit weight water.
        /// </summary>
        /// <value>
        /// The unit weight water.
        /// </value>
        public double UnitWeightWater { get; set; }

        /// <summary>
        /// Gets or sets a value indicating whether this <see cref="PlLinesCreator"/> is inwards.
        /// </summary>
        /// <value>
        ///   <c>true</c> if inwards; otherwise, <c>false</c>.
        /// </value>
        public bool Inwards { get; set; }

        /// <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 CreatePLLine(PLLineType plLineType, Location location)
        {
            GeometryPointString plLine = null;

            switch (plLineType)
            {
                case PLLineType.PL1:
                    plLine = CreatePLLine1(location);
                    break;
                case PLLineType.PL2:
                    plLine = CreatePLLine2(location);
                    break;
                case PLLineType.PL3:
                    plLine = CreatePLLine3(location);
                    break;
                case PLLineType.PL4:
                    plLine = CreatePLLine4(location);
                    break;
            }

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

        /// <summary>
        /// Create PL2 (i.e. PL-line at the penetration zone level in clay)
        /// </summary>
        /// <returns>PL-line 2 (also called Head 2)</returns>
        public GeometryPointString CreatePLLine2(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 = CreatePLLine2For1DGeometry(location);
            }

            return plLine;
        }

        /// <summary>
        /// Create PL-line 3 (= head in the Pleistocene)
        /// </summary>
        /// <returns>PL-line 3</returns>
        public GeometryPointString CreatePLLine3(Location location)
        {
            switch (location.DikeSoilScenario)
            {
                case DikeSoilScenario.ClayDikeOnClay:
                case DikeSoilScenario.SandDikeOnClay:
                    return CreatePLLine3Or4ForDikesOnClay(location.HeadInPLLine3, PLLineType.PL3, location);
                case DikeSoilScenario.ClayDikeOnSand:
                    return createPLLine3ForClayDikeOnSand(location);
            }
            return null;
        }

        /// <summary>
        /// Create PL-line 4 (= head in the in-between aquifer)
        /// </summary>
        /// <returns>PL-line 4</returns>
        public GeometryPointString CreatePLLine4(Location location)
        {
            switch (location.DikeSoilScenario)
            {
                case DikeSoilScenario.ClayDikeOnClay:
                case DikeSoilScenario.SandDikeOnClay:
                    return CreatePLLine3Or4ForDikesOnClay(location.HeadInPLLine4, PLLineType.PL4, location);
            }
            return null;
        }

        /// <summary>
        /// Determines the intersection between two line segments:
        /// - Horizontal water level
        /// - Surface line (until DikeTopAtRiver)
        /// </summary>
        /// <param name="waterLevelRiver">Water Level at River side</param>
        /// <param name="surfaceline">The surface line</param>
        /// <returns>Intersection point between horizontal water level and surface line</returns>
        public GeometryPoint DetermineIntersectionBetweenWaterLevelAndDike(double waterLevelRiver,
                                                                           SurfaceLine2 surfaceline)
        {
            GeometryPoint pointAtLevel = DetermineIntersectionBetweenTaludRiverSideAndWaterLevel(waterLevelRiver,
                                                                                                 surfaceline);
            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;
            }
        }

        /// <summary>
        /// Determines the intersection between PL2 or PL1 and the final slope (1:100) of PL3.
        /// </summary>
        /// <param name="firstPointOfObliqueLine">The first point of oblique line.</param>
        /// <param name="slopeOfObliqueLine">The slope of oblique line.</param>
        /// <param name="location">The location.</param>
        /// <returns>The intersection point.</returns>
        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 (Pl2 != null)
            {
                if (CalculateHeadPl2AtX(location, firstPointOfObliqueLine.X) < firstPointOfObliqueLine.Z)
                {
                    signSlope = -1;
                }
                else
                {
                    signSlope = 1;
                }
            }
            else
            {
                if (phreaticLine.GetZAtX(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;
        }

        /// <summary>
        /// Determines if values <paramref name="aValue1"/> and <paramref name="aValue2"/> are equal, with a tolerance of <see cref="cEpsilon"/>.
        /// </summary>
        /// <param name="aValue1">a value1.</param>
        /// <param name="aValue2">a value2.</param>
        /// <returns><c>true</c> if values <paramref name="aValue1"/> and <paramref name="aValue2"/> are equal; otherwise <c>false</c>.</returns>
        public static bool Equal(double aValue1, double aValue2)
        {
            if (Math.Abs(aValue1 - aValue2) < cEpsilon)
            {
                return true;
            }

            return false;
        }

        /// <summary>
        /// Create PL1 (= phreatic level)
        /// </summary>
        /// <returns>PL1 (= phreatic level)</returns>
        public GeometryPointString CreatePLLine1(Location location)
        {
            GeometryPoint dikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
            GeometryPoint dikeTopAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);

            // Create Phreatic line and add polder water 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;
            }

            // 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)                                                          
            CreatePhreaticLineSegmentsInsideDikeForRingtoetsWti2017(line, location);
            CreatePhreaticLineSegmentsInShoulderAndPolderForRingtoetsWti2017(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

            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);

            // 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);
        }

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

        /// <summary>
        /// Get soil profile 1D below <paramref name="xCoordinate"/> using <paramref name="soilProfile2D"/>.
        /// </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 CreatePLLine2For1DGeometry(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;
        }

        private GeometryPointString createPLLine3ForClayDikeOnSand(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 CreatePLLine3Or4ForDikesOnClay(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;
                }

                // 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 Top 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;
            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)
        {
            SoilProfile1D actualSoilProfile = location.SoilProfile2D != null
                                                  ? GetSoilProfileBelowPoint(location.SoilProfile2D, xCoor)
                                                  : location.SoilProfile1D;
            SoilLayer1D relevantSandLayer = null;
            try
            {
                relevantSandLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile, location);
            }
            catch
            {
                return Double.MaxValue;
            }

            UpliftCalculatorWTI upliftCalculator = new UpliftCalculatorWTI();
            upliftCalculator.PhreaticLevel = phreaticLine.GetZAtX(xCoor);
            upliftCalculator.SoilProfile = actualSoilProfile;
            upliftCalculator.SurfaceLevel = actualSoilProfile.TopLevel;
            upliftCalculator.TopOfLayerToBeEvaluated = relevantSandLayer.TopLevel;
            upliftCalculator.UpLiftTopLevel = relevantSandLayer.TopLevel;
            upliftCalculator.VolumicWeightOfWater = UnitWeightWater;

            return upliftCalculator.CalculateHeadOfPLLine(cUpliftFactorEquilibrium);
        }

        /// <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 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();
        }

        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.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>
        /// 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>
        /// 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.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;
                }
            }
        }

        /// <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">The phreatic line.</param>
        /// <param name="location">The location.</param>
        /// <System.Exception cref="Exception"></System.Exception>
        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 System.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);
            }
        }
    }
}