// Copyright (C) Stichting Deltares 2024. All rights reserved.
// 
// This file is part of the Dam Engine.
// 
// The Dam Engine is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero 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 Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero 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.DamEngine.Calculators.KernelWrappers.Common;
using Deltares.DamEngine.Calculators.Properties;
using Deltares.DamEngine.Calculators.Uplift;
using Deltares.DamEngine.Data.General;
using Deltares.DamEngine.Data.General.Gauges;
using Deltares.DamEngine.Data.General.PlLines;
using Deltares.DamEngine.Data.Geometry;
using Deltares.DamEngine.Data.Geotechnics;
using Deltares.DamEngine.Data.Standard;

namespace Deltares.DamEngine.Calculators.PlLinesCreator;

/// <summary>
/// Creator for Piezometric level lines
/// </summary>
public class PlLinesCreator
{
    private const double toleranceAlmostEquals = 1e-6;
    private const double cUpliftFactorEquilibrium = 1.0;
    private const double cOffsetPhreaticLineBelowSurface = 0.01;
    protected readonly Dictionary<double, SoilProfile1D> cachedSoilProfiles1D = new Dictionary<double, SoilProfile1D>();

    private double? waterLevelRiverLow;

    internal PlLine CurrentPl1Line; // is needed when calculating uplift reduction for PL3 and Pl4
    
    /// <summary>
    /// Constructor
    /// </summary>
    public PlLinesCreator()
    {
        SoilProfileType = SoilProfileType.ProfileType1D;
        PlLineOffsetBelowDikeTopAtRiver = 0.5; // Default value
        PlLineOffsetBelowDikeTopAtPolder = 1.5; // Default value
        PlLineOffsetBelowShoulderBaseInside = 0.1; // Default value
        PlLineOffsetBelowDikeToeAtPolder = 0.1; // Default value            
    }

    public SoilProfileType SoilProfileType { get; set; }
    public SoilProfile1D SoilProfile { get; set; }
    public SoilProfile2D SoilProfile2D { get; set; }
    public Soil DikeEmbankmentMaterial { get; set; }
    public SoilList SoilList { get; set; }
    public double? HeadInPlLine2 { get; set; }
    public bool IsAdjustPL3AndPL4SoNoUpliftWillOccurEnabled { get; set; } = true;
    public double WaterLevelRiverHigh { get; set; }

    public double? WaterLevelRiverLow
    {
        get
        {
            return waterLevelRiverLow;
        }
        set
        {
            waterLevelRiverLow = value;
        }
    }

    public bool IsUseLowWaterLevel { get; set; }

    public double WaterLevelPolder { get; set; }

    public double? HeadInPlLine3
    {
        get
        {
            return headInPlLine3;
        }
        set
        {
            headInPlLine3 = value;
        }
    }

    public double? HeadInPlLine4
    {
        get
        {
            return headInPlLine4;
        }
        set
        {
            headInPlLine4 = value;
        }
    }

    /// <remarks>Aggregation relationship</remarks>
    public SurfaceLine2 SurfaceLine { get; set; }

    public IList<GaugePlLine> GaugePlLines { get; set; }
    public IList<Gauge> Gauges { get; set; }
    public double GaugeMissVal { get; set; }
    public double PlLineOffsetBelowDikeTopAtRiver { get; set; }
    public double PlLineOffsetBelowDikeTopAtPolder { get; set; }
    public virtual double PlLineOffsetBelowShoulderBaseInside { get; set; }
    public virtual double PlLineOffsetBelowDikeToeAtPolder { get; set; }
    public double? PlLineOffsetBelowDikeCrestMiddle { get; set; }
    public double? PlLineOffsetFactorBelowShoulderCrest { get; set; }
    public bool? UsePlLineOffsetBelowDikeCrestMiddle { get; set; }
    public bool? UsePlLineOffsetFactorBelowShoulderCrest { get; set; }
    public double Pl3MinUplift { get; private set; }
    public double Pl3HeadAdjusted { get; private set; }
    public double Pl3LocationXMinUplift { get; private set; }
    public double Pl4MinUplift { get; private set; }
    public double Pl4HeadAdjusted { get; private set; }
    public double Pl4LocationXMinUplift { get; private set; }

    public ModelParametersForPlLines ModelParametersForPlLines { get; set; } = new ModelParametersForPlLines();

    /// <summary>
    /// Create all PlLines
    /// </summary>
    /// <returns></returns>
    public PlLines CreateAllPlLines(Location location)
    {
        var plLines = new PlLines();
        switch (ModelParametersForPlLines.PlLineCreationMethod)
        {
            case PlLineCreationMethod.ExpertKnowledgeLinearInDike:
            case PlLineCreationMethod.ExpertKnowledgeRRD:
                plLines = CreateAllPlLinesWithExpertKnowledge(location);
                break;
            case PlLineCreationMethod.GaugesWithFallbackToExpertKnowledgeRRD:
                plLines = CreateAllPlLinesWithGaugesWithFallbackToExpertKnowledgeRrd(location);
                break;
        }

        // If PL line 2 does not exist, make it equal to PL line 4
        if (!plLines.Lines[PlLineType.Pl2].Exists())
        {
            plLines.Lines[PlLineType.Pl2] = plLines.Lines[PlLineType.Pl4].Clone();
        }

        if (!plLines.Lines[PlLineType.Pl1].IsXAscending())
        {
            throw new PlLinesCreatorException("PlLine 1 not an X-ascending polyline");
        }

        return plLines;
    }

    /// <summary>
    /// Creates the pl line by expert knowledge.
    /// </summary>
    /// <param name="plLineType">Type of the pl line.</param>
    /// <param name="slopeGradient">The slope gradient.</param>
    /// <returns></returns>
    public PlLine CreatePlLineByExpertKnowledge(PlLineType plLineType, double slopeGradient)
    {
        PlLine plLine = null;

        switch (plLineType)
        {
            case PlLineType.Pl1:
                plLine = CreatePlLine1ByExpertKnowledge();
                break;
            case PlLineType.Pl2:
                plLine = CreatePlLine2ByExpertKnowledge(ModelParametersForPlLines.PenetrationLength, HeadInPlLine2);
                break;
            case PlLineType.Pl3:
                plLine = CreatePlLine3ByExpertKnowledge(DetermineHeadPl3(), ModelParametersForPlLines.DampingFactorPl3, slopeGradient);
                break;
            case PlLineType.Pl4:
                plLine = CreatePlLine4ByExpertKnowledge(DetermineHeadPl4(), ModelParametersForPlLines.DampingFactorPl4, slopeGradient);
                break;
        }

        if (plLine != null)
        {
            plLine.DeleteCoincidingPoints();
        }

        return plLine;
    }

    /// <summary>
    /// Intersection between two line segments:
    /// - Horizontal water level
    /// - Surface line (until DikeTopAtRiver)
    /// </summary>
    /// <param name="waterLevelRiver"></param>
    /// <returns></returns>
    /// 
    private PlLinePoint DetermineIntersectionBetweenWaterLevelAndDike(double waterLevelRiver)
    {
        GeometryPoint pointAtLevel = SurfaceLine.DetermineIntersectionWithLevel(waterLevelRiver);
        if (pointAtLevel != null)
        {
            return new PlLinePoint(pointAtLevel.X, pointAtLevel.Z);
        }

        return null;
    }

    private double? headInPlLine3 { get; set; }
    private double? headInPlLine4 { get; set; }

    private double WaterLevelToUse()
    {
        if (IsUseLowWaterLevel)
        {
            return waterLevelRiverLow.Value;
        }

        return WaterLevelRiverHigh;
    }

    private double DetermineHeadPl3()
    {
        double waterLevel = WaterLevelToUse();
        // If no value known for headPl3 then use water level for headPL3
        if (HeadInPlLine3 == null)
        {
            return waterLevel;
        }

        return HeadInPlLine3.Value;
    }

    private double DetermineHeadPl4()
    {
        double waterLevel = WaterLevelToUse();
        if (HeadInPlLine4 == null)
        {
            // If no value known for headPl4 then use water level (buitenwaterstand) for headPL4
            return waterLevel;
        }

        // If value is specified for headPl4 then use this value for headPL4
        return HeadInPlLine4.Value;
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="xCoordinate"></param>
    /// <returns></returns>
    private SoilProfile1D GetSoilProfileBelowPoint(double xCoordinate)
    {
        if (cachedSoilProfiles1D.ContainsKey(xCoordinate))
        {
            return cachedSoilProfiles1D[xCoordinate];
        }

        SoilProfile1D soilProfile;
        switch (SoilProfileType)
        {
            case SoilProfileType.ProfileType1D:
                soilProfile = SoilProfileHelper.DetermineForSurfaceLineCorrected1DProfileAtX(SoilProfile, SurfaceLine, xCoordinate, DikeEmbankmentMaterial);
                cachedSoilProfiles1D[xCoordinate] = soilProfile;
                return soilProfile;
            case SoilProfileType.ProfileType2D:
                soilProfile = SoilProfile2D.GetSoilProfile1D(xCoordinate);
                cachedSoilProfiles1D[xCoordinate] = soilProfile;
                return soilProfile;
            default:
                return null;
        }
    }

    /// <summary>
    /// Create PL2 (is pl line for penetration)
    /// </summary>
    /// <returns></returns>
    private PlLine CreatePlLine2ByExpertKnowledge(double penetrationLength, double? headInPlLine2)
    {
        PlLine plLine = null;

        if (penetrationLength < 0.0)
        {
            throw new PlLinesCreatorException("Negative penetration length.");
        }

        if (penetrationLength.AlmostEquals(0.0) || (headInPlLine2 == null))
        {
            // No penetration, or no Head Pl2 defined, so empty pl-line will be returned
            plLine = new PlLine();
        }
        else
        {
            switch (SoilProfileType)
            {
                case SoilProfileType.ProfileType1D:
                    plLine = CreatePlLine2ByExpertKnowledgeFor1DGeometry(penetrationLength, headInPlLine2);
                    break;
                case SoilProfileType.ProfileType2D:
                    plLine = CreatePlLine2ByExpertKnowledgeFor2DGeometry(headInPlLine2);
                    break;
            }
        }

        return plLine;
    }

    /// <summary>
    /// Create PL2 (is pl line for penetration) for 2d-geometry
    /// </summary>
    /// <param name="headInPlLine2">The head in pl line2.</param>
    /// <returns></returns>
    /// <exception cref="PlLinesCreatorException">Head PL2 not defined</exception>
    private PlLine CreatePlLine2ByExpertKnowledgeFor2DGeometry(double? headInPlLine2)
    {
        if (headInPlLine2 == null)
        {
            throw new PlLinesCreatorException("Head PL2 not defined");
        }

        var plLine = new PlLine();
        plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.First().X, headInPlLine2.Value));
        plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.Last().X, headInPlLine2.Value));
        return plLine;
    }

    /// <summary>
    /// Create PL2 (is pl line for penetration) for 1d-geometry
    /// </summary>
    /// <param name="penetrationLength"></param>
    /// <param name="headInPlLine2"></param>
    /// <returns></returns>
    private PlLine CreatePlLine2ByExpertKnowledgeFor1DGeometry(double penetrationLength, double? headInPlLine2)
    {
        if (headInPlLine2 == null)
        {
            throw new PlLinesCreatorException("Head PL2 not defined");
        }

        var plLine = new PlLine();
        if (SoilProfile != null)
        {
            SoilProfile1D relevantSoilProfile = GetRelevantSoilProfileForAquiferLayersSearch();
            IList<SoilLayer1D> aquiferLayers = relevantSoilProfile.GetAquiferLayers();
            if (aquiferLayers.Count == 0)
            {
                throw new PlLinesCreatorException("Soil profile contains no aquifer layers.");
            }

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

                if (infiltrationLayers.Count > 0)
                {
                    double separationLevel = relevantSoilProfile.BottomAquiferLayer.TopLevel + penetrationLength;

                    if (separationLevel <= infiltrationLayers.First().TopLevel)
                    {
                        plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.First().X, headInPlLine2.Value));
                        plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.Last().X, headInPlLine2.Value));
                    }
                }
            }
        }

        return plLine;
    }

    /// <summary>
    /// Create PL3
    /// </summary>
    /// <returns></returns>
    private PlLine CreatePlLine3ByExpertKnowledge(double headValue, double dampingFactor, double slopeGradient)
    {
        return CreatePlLine3Or4ByExpertKnowledge(headValue, dampingFactor, PlLineType.Pl3, slopeGradient);
    }

    /// <summary>
    /// Create PL4
    /// </summary>
    /// <returns></returns>
    private PlLine CreatePlLine4ByExpertKnowledge(double headValue, double dampingFactor, double slopeGradient)
    {
        return CreatePlLine3Or4ByExpertKnowledge(headValue, dampingFactor, PlLineType.Pl4, slopeGradient);
    }

    private PlLine CreatePlLine3Or4ByExpertKnowledge(double headValue, double dampingFactor, PlLineType plLineType, double slopeGradient)
    {
        var plLine = new PlLine();

        if (dampingFactor < 0.0)
        {
            throw new PlLinesCreatorException("Damping factor < 0.0");
        }

        // Soil profile 1D below toe of dike at polder is used to check if there is really a relevant aquifer.
        SoilProfile1D actualSoilProfile = GetRelevantSoilProfileForAquiferLayersSearch();
        SoilLayer1D relevantAquiferLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile);
        if (relevantAquiferLayer != null)
        {
            double referenceLevel = HeadInPlLine2 ?? WaterLevelPolder;
            double headAtPolderDikeToe = headValue - Math.Max(0, dampingFactor * (headValue - referenceLevel));
            plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.First().X, headValue));
            plLine.Points.Add(new PlLinePoint(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X, headValue));

            if (IsSurfaceLineIntersectedByAquiferAtPolder(plLineType, out GeometryPoint intersectionAquifer))
            {
                ReducePlLineToPl1(intersectionAquifer, plLine);
            }
            else
            {
                plLine.Points.Add(new PlLinePoint(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, headAtPolderDikeToe));

                // Now continue PlLine to the end with a slope of SlopeGradient
                AddTailOfPl3OrPl4WithSlopeGradient(slopeGradient, plLine);

                if (IsAdjustPL3AndPL4SoNoUpliftWillOccurEnabled)
                {
                    AdjustLineAccordingToTrwUplift(plLine, plLineType, slopeGradient);
                }
            }

            EnsureDescendingLine(plLine);

            RemoveRedundantPoints(plLine);
        }

        return plLine;
    }

    /// <summary>
    /// Refer to Word document "20250116 Schematisatie freatische lijn en waterspanningen waarbij teensloot een aquifer insnijdt"
    /// in doc/Work folder for details on the reduction of the PL-line.
    /// </summary>
    /// <param name="intersectionAquifer">The intersection point of the surface line (often ditch) with the aquifer.</param>
    /// <param name="plLine">The PL-line.</param>
    internal void ReducePlLineToPl1(GeometryPoint intersectionAquifer, PlLine plLine)
    {
        if (IsDitchIntersectedByPl1(intersectionAquifer, out GeometryPoint intersectDitchPl1))
        {
            plLine.Points.Add(new PlLinePoint(intersectDitchPl1.X, intersectDitchPl1.Z));
            plLine.Points.Add(new PlLinePoint(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X, intersectDitchPl1.Z));
        }
        else
        {
            double pl1Level = CurrentPl1Line.ZFromX(intersectionAquifer.X);
            plLine.Points.Add(new PlLinePoint(intersectionAquifer.X, pl1Level));
            plLine.Points.Add(new PlLinePoint(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X, pl1Level));
        }
    }

    /// <summary>
    ///  Continue PlLine to the end with a slope of slope gradient
    ///  If PlLine is lower than polder level, continue with polder level
    /// </summary>
    /// <param name="slopeGradient">The slope gradient.</param>
    /// <param name="plLine">The pl line.</param>
    private void AddTailOfPl3OrPl4WithSlopeGradient(double slopeGradient, PlLine plLine)
    {
        if (plLine.Points.Last().Z <= WaterLevelPolder)
        {
            // the PL line is already at WaterLevelPolder, so countinue it horizontally
            plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.Last().X, WaterLevelPolder));
        }
        else
        {
            double lengthFromLastPlPointToEnd = Math.Abs(SurfaceLine.Geometry.Points.Last().X - plLine.Points.Last().X);
            double headAtLastPlPoint = plLine.Points.Last().Z;
            double headAtEnd = plLine.Points.Last().Z - slopeGradient * lengthFromLastPlPointToEnd;

            var waterLevelPolderLine =
                new Line<Point2D>(new Point2D(SurfaceLine.Geometry.Points.First().X, WaterLevelPolder),
                                  new Point2D(SurfaceLine.Geometry.Points.Last().X, WaterLevelPolder));
            var slopeLine =
                new Line<Point2D>(
                    new Point2D(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, headAtLastPlPoint),
                    new Point2D(SurfaceLine.Geometry.Points.Last().X, headAtEnd));

            if (waterLevelPolderLine.IntersectsZ(slopeLine, out Point2D intersectionPoint))
            {
                plLine.Points.Add(new PlLinePoint(intersectionPoint.X, WaterLevelPolder));
                plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.Last().X, WaterLevelPolder));
            }
            else
            {
                plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.Last().X, headAtEnd));
            }
        }
    }

    private SoilLayer1D GetRelevantAquiferLayer(PlLineType type, SoilProfile1D actualSoilProfile)
    {
        SoilLayer1D relevantAquiferLayer;
        switch (type)
        {
            case PlLineType.Pl3:
                relevantAquiferLayer = actualSoilProfile.BottomAquiferLayer; // Pleistocene
                break;
            case PlLineType.Pl4:
                relevantAquiferLayer = actualSoilProfile.InBetweenAquiferLayer; // Intermediate sand layer 
                break;
            default:
                throw new PlLinesCreatorException($"Invalid PL line type:{type} for creation of PL Line 3 or 4");
        }

        return relevantAquiferLayer;
    }

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

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

            plPointPrevious = plPoint;
        }
    }

    /// <summary>
    /// Clear outputvalues for PL3 or PL4
    /// </summary>
    /// <param name="plLineType"></param>
    private void ClearOutputValuesForPl3_4(PlLineType plLineType)
    {
        switch (plLineType)
        {
            case PlLineType.Pl3:
                Pl3MinUplift = Double.MaxValue;
                Pl3HeadAdjusted = Double.MaxValue;
                Pl3LocationXMinUplift = Double.MaxValue;
                break;
            case PlLineType.Pl4:
                Pl4MinUplift = Double.MaxValue;
                Pl4HeadAdjusted = Double.MaxValue;
                Pl4LocationXMinUplift = Double.MaxValue;
                break;
        }
    }

    /// <summary>
    /// Determine outputvalues for PL3 or PL4 for minimal upliftfactor
    /// </summary>
    /// <param name="plLineType">Type of the pl line.</param>
    /// <param name="xCoordinate">The x coordinate.</param>
    /// <param name="upliftFactor">The uplift factor.</param>
    /// <param name="headValue">The head value.</param>
    private void UpdateOutputValuesForPl3_4(PlLineType plLineType, double xCoordinate, double upliftFactor, double headValue)
    {
        switch (plLineType)
        {
            case PlLineType.Pl3:
                if (upliftFactor < Pl3MinUplift)
                {
                    Pl3MinUplift = upliftFactor;
                    Pl3HeadAdjusted = headValue;
                    Pl3LocationXMinUplift = xCoordinate;
                }

                break;
            case PlLineType.Pl4:
                if (upliftFactor < Pl4MinUplift)
                {
                    Pl4MinUplift = upliftFactor;
                    Pl4HeadAdjusted = headValue;
                    Pl4LocationXMinUplift = xCoordinate;
                }

                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 FindIndexOfPointInPlLineWithXCoordinate(PlLine plLine, double x)
    {
        for (int pointIndex = plLine.Points.Count - 1; pointIndex > 0; pointIndex--)
        {
            PlLinePoint plPoint = plLine.Points[pointIndex];
            if (plPoint.X.AlmostEquals(x, GeometryPoint.Precision))
            {
                return pointIndex;
            }
        }

        return -1;
    }

    /// <summary>
    /// Removes the index of all points of pl line after the start index.
    /// </summary>
    /// <param name="plLine">The pl line.</param>
    /// <param name="startIndex">The start index (this point will not be removed).</param>
    private void RemoveAllPointsOfPlLineAfterIndex(PlLine plLine, int startIndex)
    {
        for (int pointIndex = plLine.Points.Count - 1; pointIndex > startIndex; pointIndex--)
        {
            plLine.Points.RemoveAt(pointIndex);
        }
    }

    /// <summary>
    /// Determines whether a x coordinate is in ditch (excluding first and last point of ditch).
    /// </summary>
    /// <param name="xCoordinate">The x coordinate.</param>
    /// <returns></returns>
    private bool IsXCoordinateInDitch(double xCoordinate)
    {
        if (SurfaceLine.HasDitch())
        {
            return (xCoordinate > SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X &&
                    xCoordinate < SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X);
        }

        // if ditch does not exist the x coordinate is clearly not in the ditch
        return false;
    }

    /// <summary>
    /// Determines the thickness of aquitard at characteristic point.
    /// This is the distance between the surface level and the top of the first aquifer
    /// </summary>
    /// <param name="characteristicPointType">Type of the characteristic point.</param>
    /// <returns></returns>
    private double? DetermineThicknessAquitardAtCharacteristicPoint(CharacteristicPointType characteristicPointType)
    {
        if (!SurfaceLine.HasAnnotation(characteristicPointType))
        {
            return null;
        }

        GeometryPoint characteristicGeometryPoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(characteristicPointType);
        SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(characteristicGeometryPoint.X);
        double bottomLevel = actualSoilProfile.GetHighestAquifer().TopLevel;
        double topLevel = characteristicGeometryPoint.Z;
        return topLevel - bottomLevel;
    }

    /// <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 slope gradient until polder level.
    /// 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 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="slopeGradient"> </param>
    private void AdjustLineAccordingToTrwUplift(PlLine plLine, PlLineType plLineType, double slopeGradient)
    {
        const double cTolerancePoint = 0.0001;
        ClearOutputValuesForPl3_4(plLineType);

        var upliftCalculator = new UpliftCalculator
        {
            IsUseOvenDryUnitWeight = false,
            UnitWeightSoilEmbankment = (DikeEmbankmentMaterial == null) ? null : DikeEmbankmentMaterial.AbovePhreaticLevel
        };
        GeometryPoint startSurfacePoint = SurfaceLine.GetDikeToeInward();

        int indexOfFixedPlPoint = FindIndexOfPointInPlLineWithXCoordinate(plLine, SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X);
        if (indexOfFixedPlPoint < 0)
        {
            throw new PlLinesCreatorException("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 = from GeometryPoint point in SurfaceLine.Geometry.Points
                                                               where point.X >= startSurfacePoint.X
                                                               orderby point.X
                                                               select point;
        // Adjustment will only be applied if the value to adjust to is smaller than the previous adjusted value
        // (to avoid that the PL3/PL4 will be adjusted to the polder side of the ditch i.s.o. the dike side of the ditch).
        // So we remember which was the last adjusted value in lastAdjustedHeadOfPlLine.
        var lastAdjustedHeadOfPlLine = Double.MaxValue;

        bool isRemoveAllPlPointsBackToDikeToeAtRiver;
        bool isSkipAdjustmentInDitch;
        DetermineHowToActDueToDitch(out isSkipAdjustmentInDitch, out isRemoveAllPlPointsBackToDikeToeAtRiver);

        foreach (GeometryPoint surfacePoint in relevantSurfacePointsList)
        {
            if (IsXCoordinateInDitch(surfacePoint.X) && isSkipAdjustmentInDitch)
            {
                continue;
            }

            ConfigureUpliftCalculator(upliftCalculator, surfacePoint);
            SoilLayer1D relevantSandLayer = GetRelevantAquiferLayer(plLineType, upliftCalculator.SoilProfile);
            // When PlLineType is PL4 it is possible that no relevant aquifer layer is found, then the adjustment is not necessary
            if (relevantSandLayer != null)
            {
                double aquiferTopLayer = Math.Min(surfacePoint.Z, relevantSandLayer.TopLevel);
                upliftCalculator.TopOfLayerToBeEvaluated = aquiferTopLayer;
                double headOfPlLine = plLine.ZFromX(surfacePoint.X);
                double upliftFactor = upliftCalculator.CalculateUpliftFactor(headOfPlLine);
                if (upliftFactor <= cUpliftFactorEquilibrium)
                {
                    // Adjust headOfPlLine so there is equilibrium
                    bool hasNoCoverLayer = surfacePoint.Z <= aquiferTopLayer;
                    if (hasNoCoverLayer)
                    {
                        headOfPlLine = WaterLevelPolder;
                    }
                    else
                    {
                        headOfPlLine = Math.Max(upliftCalculator.CalculateHeadOfPlLine(cUpliftFactorEquilibrium),
                                                WaterLevelPolder);
                    }

                    if (headOfPlLine < lastAdjustedHeadOfPlLine)
                    {
                        var plPoint = new PlLinePoint(surfacePoint.X, headOfPlLine);
                        // Remove all points of PlLine after fixed point
                        RemoveAllPointsOfPlLineAfterIndex(plLine, indexOfFixedPlPoint);
                        plLine.Points.Add(plPoint);
                        if (!isRemoveAllPlPointsBackToDikeToeAtRiver)
                        {
                            // To make sure that not all points of the PL-line back to the toe of the dike at riverside are to be removed
                            indexOfFixedPlPoint = plLine.Points.Count - 1;
                        }

                        AddTailOfPl3OrPl4WithSlopeGradient(slopeGradient, plLine);
                        lastAdjustedHeadOfPlLine = headOfPlLine;
                    }
                }

                UpdateOutputValuesForPl3_4(plLineType, surfacePoint.X, upliftFactor, headOfPlLine);
            }
        }

        // Recheck on Uplift in case points are removed 
        if (isRemoveAllPlPointsBackToDikeToeAtRiver && SurfaceLine.HasDitch())
        {
            GeometryPoint startPoint = SurfaceLine.GetDikeToeInward();
            GeometryPoint endPoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide);
            IEnumerable<GeometryPoint> recheckSurfacePointsList = from GeometryPoint point in SurfaceLine.Geometry.Points
                                                                  where point.X >= startPoint.X && point.X <= endPoint.X
                                                                  orderby point.X
                                                                  select point;

            foreach (GeometryPoint surfacePoint in recheckSurfacePointsList)
            {
                ConfigureUpliftCalculator(upliftCalculator, surfacePoint);
                SoilLayer1D relevantSandLayer = GetRelevantAquiferLayer(plLineType, upliftCalculator.SoilProfile);
                // When PlLineType is PL4 it is possible that no relevant aquifer layer is found, then the adjustment is not necessary
                if (relevantSandLayer != null)
                {
                    double aquiferTopLayer = Math.Min(surfacePoint.Z, relevantSandLayer.TopLevel);
                    upliftCalculator.TopOfLayerToBeEvaluated = aquiferTopLayer;
                    double headOfPlLine = plLine.ZFromX(surfacePoint.X);
                    double upliftFactor = upliftCalculator.CalculateUpliftFactor(headOfPlLine);
                    if (upliftFactor <= cUpliftFactorEquilibrium)
                    {
                        // Adjust headOfPlLine so there is equilibrium
                        headOfPlLine = Math.Max(upliftCalculator.CalculateHeadOfPlLine(cUpliftFactorEquilibrium), WaterLevelPolder);
                        double currentHead = plLine.ZFromX(surfacePoint.X);
                        if (headOfPlLine < currentHead)
                        {
                            PlLinePoint plPoint = plLine.EnsurePointAtX(surfacePoint.X, cTolerancePoint);
                            plPoint.Z = headOfPlLine;
                        }
                    }

                    UpdateOutputValuesForPl3_4(plLineType, surfacePoint.X, upliftFactor, headOfPlLine);
                }
            }
        }
    }

    /// <summary>
    /// Configures the uplift calculator.
    /// </summary>
    /// <param name="upliftCalculator">The uplift calculator.</param>
    /// <param name="surfacePoint">The surface point.</param>
    private void ConfigureUpliftCalculator(UpliftCalculator upliftCalculator, GeometryPoint surfacePoint)
    {
        upliftCalculator.SurfaceLevel = surfacePoint.Z;
        if (CurrentPl1Line != null)
        {
            upliftCalculator.PhreaticLevel = CurrentPl1Line.ZFromX(surfacePoint.X);
            // set phreatic level to calculate uplift factor
        }
        else
        {
            upliftCalculator.PhreaticLevel = WaterLevelPolder;
            // if not PL1 created then assume phreatic level is same as waterlevel polder
        }

        SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X);
        AdaptSoilProfileForSurfaceLevelAccuracy(actualSoilProfile, upliftCalculator.SurfaceLevel);
        upliftCalculator.SoilProfile = actualSoilProfile;
    }

    /// <summary>
    /// Determines how to act due to ditch.
    /// </summary>
    /// <param name="isSkipAdjustmentInDitch">if set to <c>true</c> [is skip adjustment in ditch].</param>
    /// <param name="isRemoveAllPlPointsBackToDikeToeAtRiver">if set to <c>true</c> [is remove all pl points back to dike toe at river].</param>
    private void DetermineHowToActDueToDitch(out bool isSkipAdjustmentInDitch, out bool isRemoveAllPlPointsBackToDikeToeAtRiver)
    {
        if (SurfaceLine.HasDitch())
        {
            double thicknessAquitardAtTopDitch =
                DetermineThicknessAquitardAtCharacteristicPoint(CharacteristicPointType.DitchPolderSide).Value;
            double thicknessAquitardAtBottomDitch =
                DetermineThicknessAquitardAtCharacteristicPoint(CharacteristicPointType.BottomDitchPolderSide).Value;
            double widthDitchAtTop = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X -
                                     SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X;
            double widthDitchAtTBottom =
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchPolderSide).X -
                SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchDikeSide).X;
            if ((thicknessAquitardAtTopDitch > widthDitchAtTop) &&
                (thicknessAquitardAtBottomDitch > widthDitchAtTBottom))
            {
                isRemoveAllPlPointsBackToDikeToeAtRiver = false;
                isSkipAdjustmentInDitch = true;
            }
            else
            {
                isRemoveAllPlPointsBackToDikeToeAtRiver = true;
                isSkipAdjustmentInDitch = false;
            }
        }
        else
        {
            isRemoveAllPlPointsBackToDikeToeAtRiver = false;
            isSkipAdjustmentInDitch = false;
        }
    }

    /// <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 void AdaptSoilProfileForSurfaceLevelAccuracy(SoilProfile1D actualSoilProfile, double surfaceLevel)
    {
        if (surfaceLevel > actualSoilProfile.TopLevel)
        {
            actualSoilProfile.Layers[0].TopLevel = surfaceLevel;
        }
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="plLine"></param>
    /// <param name="genericLine"></param>
    void CopyPointsInPlLine(ref PlLine plLine, GeometryPointString genericLine)
    {
        plLine.Points.Clear();
        foreach (GeometryPoint point in genericLine.Points)
        {
            plLine.Points.Add(new PlLinePoint(point.X, point.Z));
        }
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="location"></param>
    /// <returns></returns>
    private PlLines CreateAllPlLinesWithExpertKnowledge(Location location)
    {
        ValidateSoilProfileForPlLinesCreation();

        var plLines = new PlLines();
        foreach (PlLineType plLineType in Enum.GetValues(typeof(PlLineType)))
        {
            if (location != null)
            {
                bool isPl1LineDefinedForLocation = location.LocalXzpl1Line != null &&
                                                   location.LocalXzpl1Line.Points.Count > 1;
                if ((plLineType == PlLineType.Pl1) && isPl1LineDefinedForLocation)
                {
                    PlLine plLine = plLines.Lines[plLineType];
                    CopyPointsInPlLine(ref plLine, location.LocalXzpl1Line);
                }
                else
                {
                    plLines.Lines[plLineType] = CreatePlLineByExpertKnowledge(plLineType,
                                                                              location.SlopeDampingPiezometricHeightPolderSide);
                }
            }

            // currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
            if (plLineType == PlLineType.Pl1)
            {
                CurrentPl1Line = plLines.Lines[plLineType];
            }
        }

        return plLines;
    }

    private PlLines CreateAllPlLinesWithGaugesWithFallbackToExpertKnowledgeRrd(Location location)
    {
        var plLines = new PlLines();

        foreach (PlLineType plLineType in Enum.GetValues(typeof(PlLineType)))
        {
            GaugePlLine gaugePlLine = (GaugePlLines != null) ? (GaugePlLines.FirstOrDefault(x => x.PlLineType == plLineType)) : null;
            if (gaugePlLine != null && location != null)
            {
                Boolean isUseWaterLevel = ((plLineType == PlLineType.Pl1) || (plLineType == PlLineType.Pl2));
                plLines.Lines[plLineType] = CreatePlLineFromGauges(gaugePlLine, Gauges, location, isUseWaterLevel);
            }
            else if (location != null)
            {
                plLines.Lines[plLineType] = CreatePlLineByExpertKnowledge(plLineType, location.SlopeDampingPiezometricHeightPolderSide);
            }

            // currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
            if (plLineType == PlLineType.Pl1)
            {
                CurrentPl1Line = plLines.Lines[plLineType];
            }
        }

        return plLines;
    }

    /// <summary>
    /// Creates the pl line from gauges.
    /// </summary>
    /// <param name="gaugePlLine">The gauge pl line.</param>
    /// <param name="gauges">The gauges.</param>
    /// <param name="location">The location.</param>
    /// <param name="useWaterLevel">if set to <c>true</c> [use water level].</param>
    /// <returns></returns>
    /// <exception cref="PlLinesCreatorException">
    /// </exception>
    private PlLine CreatePlLineFromGauges(GaugePlLine gaugePlLine, IEnumerable<Gauge> gauges, Location location, Boolean useWaterLevel)
    {
        const double toleranceAlmostEqual = 1e-09;
        var plLine = new PlLine();

        double? gaugeWaterLevelRiver = null;
        double? leftMostXAtRiverLevel = null;

        var pointIndex = 0;
        foreach (GaugePlLinePoint gaugePlLinePoint in gaugePlLine.Points)
        {
            double? localX = gaugePlLinePoint.X;
            double? localZ = gaugePlLinePoint.Z;

            if (gauges != null)
            {
                if (gaugePlLinePoint.GaugeIDX != null && gaugePlLinePoint.GaugeIDX != "")
                {
                    Gauge gauge = (gauges.Where(x => x.Name == gaugePlLinePoint.GaugeIDX && x.Location == location)).FirstOrDefault();
                    if (gauge != null)
                    {
                        localX = gauge.LocalX;
                    }
                    else
                    {
                        throw new PlLinesCreatorException($"Gauge PL line {gaugePlLine.PlLineType.ToString()} refers to an unknown gauge named '{gaugePlLinePoint.GaugeIDX}' at X coordinate #{pointIndex}.");
                    }
                }

                if (gaugePlLinePoint.GaugeIDZ != null && gaugePlLinePoint.GaugeIDZ != "")
                {
                    Gauge gauge = (gauges.Where(x => x.Name == gaugePlLinePoint.GaugeIDZ && x.Location == location)).FirstOrDefault();
                    if (gauge != null)
                    {
                        if (!gauge.Value.HasValue || ((double) gauge.Value).IsNearEqual(GaugeMissVal, toleranceAlmostEqual))
                        {
                            throw new PlLinesCreatorException($"Value of gauge {gauge.Name} at location {location.Name} in gauge PL line of type {gaugePlLine.PlLineType.ToString()} is undefined.");
                        }

                        localZ = gauge.Value;
                    }
                    else
                    {
                        throw new PlLinesCreatorException($"Gauge PL line {gaugePlLine.PlLineType.ToString()} refers to an unknown gauge named '{gaugePlLinePoint.GaugeIDZ}' at Z coordinate #{pointIndex}.");
                    }
                }
            }

            if (!localX.HasValue)
            {
                throw new PlLinesCreatorException($"Gauge PL line {gaugePlLine.PlLineType.ToString()}'s X coordinate #{pointIndex} is undefined.");
            }

            if (!localZ.HasValue)
            {
                throw new PlLinesCreatorException($"Gauge PL line {gaugePlLine.PlLineType.ToString()}'s value #{pointIndex} is undefined.");
            }

            if (!leftMostXAtRiverLevel.HasValue || localX > leftMostXAtRiverLevel)
            {
                plLine.Points.Add(new PlLinePoint(localX.Value, localZ.Value));
                if (useWaterLevel)
                {
                    // Have to account for waterlevel
                    if (!gaugeWaterLevelRiver.HasValue)
                    {
                        // Use first gauge as waterlevel
                        gaugeWaterLevelRiver = localZ.Value;
                        IList<double> intersectionsAtX = SurfaceLine.Geometry.IntersectionsXAtZ(gaugeWaterLevelRiver.Value);
                        if (intersectionsAtX.Count() > 0)
                        {
                            // this is where the waterlevel intersects with the dike
                            leftMostXAtRiverLevel = intersectionsAtX.First();
                            plLine.Points.Add(new PlLinePoint(leftMostXAtRiverLevel.Value, gaugeWaterLevelRiver.Value));
                        }
                        else
                        {
                            // No intersection with dike; polder is flooded
                            leftMostXAtRiverLevel = SurfaceLine.Geometry.Points.Last().X;
                            plLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points.Last().X, gaugeWaterLevelRiver.Value));
                        }
                    }
                }
                else
                {
                    // In this case no intersection for this PlLine with the dike will be considered
                    leftMostXAtRiverLevel = SurfaceLine.Geometry.Points.First().X;
                }
            }

            pointIndex++;
        }

        if (plLine.Points.Any())
        {
            plLine.Points.First().X = SurfaceLine.Geometry.Points.First().X;
            plLine.Points.Last().X = SurfaceLine.Geometry.Points.Last().X;
        }

        return plLine;
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="phreaticLine"></param>
    /// <param name="lowWaterLevel"></param>
    /// <param name="highWaterLevel"></param>
    private void CreatePhreaticLineSegmentsInsideDikeForLowRiverLevel(PlLine phreaticLine, double lowWaterLevel, double highWaterLevel)
    {
        PlLinePoint intersectionLowWaterLevelWithDike = DetermineIntersectionBetweenWaterLevelAndDike(lowWaterLevel);
        GeometryPoint pointDikeToeAtRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
        switch (ModelParametersForPlLines.PlLineCreationMethod)
        {
            case PlLineCreationMethod.ExpertKnowledgeLinearInDike:
                const double cPlLineOffsetBelowSurface = 0.1;
                PlLinePoint intersectionHighWaterLevelWithDike = DetermineIntersectionBetweenWaterLevelAndDike(highWaterLevel);
                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 PlLinePoint(pointDikeToeAtRiver.X, pointDikeToeAtRiver.Z);
                }

                // At this point both intersectionHighWaterLevelWithDike and intersectionLowWaterLevelWithDike are <> null
                // or else an exception would have been thrown
                var pointBelowHighLevel = new PlLinePoint(intersectionHighWaterLevelWithDike.X,
                                                          intersectionHighWaterLevelWithDike.Z - cPlLineOffsetBelowSurface);

                //Add points below surface of dike talud riverside
                foreach (GeometryPoint point in SurfaceLine.Geometry.Points.Where(
                             point => point.X > intersectionLowWaterLevelWithDike.X &&
                                      point.X < intersectionHighWaterLevelWithDike.X))
                {
                    phreaticLine.Points.Add(new PlLinePoint(point.X, 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 (GeometryPoint point in SurfaceLine.Geometry.Points.Where(
                             point => point.X > intersectionLowWaterLevelWithDike.X &&
                                      point.X <= pointDikeToeAtRiver.X))
                {
                    //if (!surfaceLine.IsNonWaterRetainingObjectPoint(point))
                    //{
                    phreaticLine.Points.Add(new PlLinePoint(point.X, point.Z - cPlLineOffsetBelowSurface));
                    //}
                }

                // Add points below crest of dike
                double offsetDikeTopAtRiver = WaterLevelRiverHigh - PlLineOffsetBelowDikeTopAtRiver;
                double offsetDikeTopAtPolder = WaterLevelRiverHigh - PlLineOffsetBelowDikeTopAtPolder;
                GeometryPoint pointDikeTopAtRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
                GeometryPoint pointDikeTopAtPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
                if (pointDikeTopAtRiver != null)
                {
                    phreaticLine.Points.Add(new PlLinePoint(pointDikeTopAtRiver.X, offsetDikeTopAtRiver));
                }

                if (pointDikeTopAtRiver != null && pointDikeTopAtPolder != null)
                {
                    CreateOptionalPointAtDikeCrestMiddle(phreaticLine, pointDikeTopAtRiver, pointDikeTopAtPolder);
                }

                if (pointDikeTopAtPolder != null)
                {
                    phreaticLine.Points.Add(new PlLinePoint(pointDikeTopAtPolder.X, offsetDikeTopAtPolder));
                }

                break;
        }
    }

    /// <summary>
    /// Creates the optional point at dike crest middle.
    /// </summary>
    /// <param name="phreaticLine">The phreatic line.</param>
    /// <param name="pointDikeTopAtRiver">The point dike top at river.</param>
    /// <param name="pointDikeTopAtPolder">The point dike top at polder.</param>
    private void CreateOptionalPointAtDikeCrestMiddle(PlLine phreaticLine, GeometryPoint pointDikeTopAtRiver,
                                                      GeometryPoint pointDikeTopAtPolder)
    {
        if (UsePlLineOffsetBelowDikeCrestMiddle.HasValue && UsePlLineOffsetBelowDikeCrestMiddle.Value && PlLineOffsetBelowDikeCrestMiddle != null)
        {
            double middleDikeCrestX = (pointDikeTopAtRiver.X + pointDikeTopAtPolder.X) * 0.5;
            double middleDikeCrestZ = WaterLevelRiverHigh - PlLineOffsetBelowDikeCrestMiddle.Value;
            // Check whether middleDikeCrestZ is above the surface line. If so, than use the value at the surfaceline instead.
            List<double> allZ = SurfaceLine.Geometry.GetAllZatXForLine(middleDikeCrestX);
            if (allZ.Count > 0)
            {
                double z = allZ.FirstOrDefault();
                if (middleDikeCrestZ > z)
                {
                    middleDikeCrestZ = z;
                }
            }

            phreaticLine.Points.Add(new PlLinePoint(middleDikeCrestX, middleDikeCrestZ));
        }
    }

    /// <summary>
    /// Create the phreatic line segments inside the dike
    /// </summary>
    /// <param name="phreaticLine"></param>
    private void CreatePhreaticLineSegmentsInsideDikeForHighRiverLevel(PlLine phreaticLine)
    {
        double offsetDikeTopAtRiver = WaterLevelRiverHigh - PlLineOffsetBelowDikeTopAtRiver;
        double offsetDikeTopAtPolder = WaterLevelRiverHigh - PlLineOffsetBelowDikeTopAtPolder;

        // points created here are in dike
        switch (ModelParametersForPlLines.PlLineCreationMethod)
        {
            case PlLineCreationMethod.ExpertKnowledgeLinearInDike:
                // No extra points below top of dike at river and top of dike at polder
                break;
            case PlLineCreationMethod.ExpertKnowledgeRRD:
                GeometryPoint pointDikeTopAtRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
                if (pointDikeTopAtRiver != null)
                {
                    phreaticLine.Points.Add(new PlLinePoint(pointDikeTopAtRiver.X, offsetDikeTopAtRiver));
                }

                GeometryPoint pointDikeTopAtPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
                if (pointDikeTopAtRiver != null && pointDikeTopAtPolder != null)
                {
                    CreateOptionalPointAtDikeCrestMiddle(phreaticLine, pointDikeTopAtRiver, pointDikeTopAtPolder);
                }

                if (pointDikeTopAtPolder != null)
                {
                    phreaticLine.Points.Add(new PlLinePoint(pointDikeTopAtPolder.X, offsetDikeTopAtPolder));
                }

                break;
        }
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="phreaticLine"></param>
    private void CreatePhreaticLineSegmentsInShoulderAndPolder(PlLine phreaticLine)
    {
        // intersectionPolderLevelWithDike is always determined between dike toe polder side and dike top riverside. Any intersection(s) at polder side
        // beyond dike toe are ignored.
        PlLinePoint intersectionPolderLevelWithDike = DetermineIntersectionBetweenPolderLevelAndDike(WaterLevelPolder);
        GeometryPoint dikeToeAtPolderPoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);

        double maxXCoordinateSurface = SurfaceLine.Geometry.Points.Max(x => x.X);
        if (ModelParametersForPlLines.PlLineCreationMethod != PlLineCreationMethod.ExpertKnowledgeLinearInDike)
        {
            // Points created below are points starting at shoulder point to the right
            GeometryPoint shoulderBaseInsidePoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside);
            if (shoulderBaseInsidePoint != null)
            {
                double zLevel = Math.Min(phreaticLine.Points.Last().Z,
                                         shoulderBaseInsidePoint.Z -
                                         Math.Max(cOffsetPhreaticLineBelowSurface, PlLineOffsetBelowShoulderBaseInside));
                zLevel = Math.Max(zLevel, WaterLevelPolder);
                // Add point if it lies left of intersection of polderlevel with dike)
                if ((intersectionPolderLevelWithDike == null) ||
                    (intersectionPolderLevelWithDike.X > shoulderBaseInsidePoint.X))
                {
                    phreaticLine.Points.Add(new PlLinePoint(shoulderBaseInsidePoint.X, zLevel));
                }

                if (UsePlLineOffsetFactorBelowShoulderCrest.HasValue && UsePlLineOffsetFactorBelowShoulderCrest.Value &&
                    PlLineOffsetFactorBelowShoulderCrest != null && dikeToeAtPolderPoint != null)
                {
                    GeometryPoint shoulderTopInsidePoint =
                        SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderTopInside);
                    if (shoulderTopInsidePoint != null)
                    {
                        zLevel = dikeToeAtPolderPoint.Z + (PlLineOffsetFactorBelowShoulderCrest.Value *
                                                           (shoulderTopInsidePoint.Z - dikeToeAtPolderPoint.Z));
                        zLevel = Math.Min(zLevel, shoulderTopInsidePoint.Z - cOffsetPhreaticLineBelowSurface);
                        zLevel = Math.Min(zLevel, phreaticLine.Points.Last().Z);
                        zLevel = Math.Max(zLevel, WaterLevelPolder);
                        phreaticLine.Points.Add(new PlLinePoint(shoulderTopInsidePoint.X, zLevel));
                    }
                }
            }
        }

        GeometryPoint ditchDikeSidePoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
        if (dikeToeAtPolderPoint != null)
        {
            double zLevel = Math.Min(phreaticLine.Points.Last().Z, dikeToeAtPolderPoint.Z - Math.Max(cOffsetPhreaticLineBelowSurface, PlLineOffsetBelowDikeToeAtPolder));
            if (ditchDikeSidePoint != null)
            {
                if (ditchDikeSidePoint.LocationEquals(dikeToeAtPolderPoint))
                {
                    // If DikeToeAtPolder is same as DitchDikeSide pl1 should always go to polderlevel at this point
                    zLevel = WaterLevelPolder;
                }
            }

            zLevel = Math.Max(zLevel, WaterLevelPolder);
            // Add point if it lies left of intersection of polderlevel with dike
            // Note Bka: this can only happen when intersectionPolderLevelWithDike = nul as dikeToeAtPolderPoint can NEVER be
            // left of intersectionPolderLevelWithDike (only when dike is totally wrong!)
            if ((intersectionPolderLevelWithDike == null) || (intersectionPolderLevelWithDike.X > dikeToeAtPolderPoint.X))
            {
                phreaticLine.Points.Add(new PlLinePoint(dikeToeAtPolderPoint.X, zLevel));
            }
        }

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

        bool isDitchPresent = (ditchDikeSidePoint != null);

        // if there is a ditch, then adjust it.
        if (isDitchPresent)
        {
            int surfacePointIndex = SurfaceLine.Geometry.Points.IndexOf(ditchDikeSidePoint);
            AdjustForDitchAtPolderSide(phreaticLine, surfacePointIndex);
        }
        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 (PlLineOffsetBelowDikeToeAtPolder > 0)
            {
                if (phreaticLine.Points[phreaticLine.Points.Count - 1].Z > WaterLevelPolder)
                {
                    AddPhreaticLineAlongSurfaceLevel(phreaticLine);
                }
            }
            else
            {
                if (dikeToeAtPolderPoint != null && phreaticLine.Points.Last().X < dikeToeAtPolderPoint.X)
                {
                    var newPoint = new PlLinePoint(dikeToeAtPolderPoint.X, WaterLevelPolder);
                    phreaticLine.Points.Add(newPoint);
                }
                else
                {
                    phreaticLine.Points.Last().Z = WaterLevelPolder;
                }
            }
        }

        //Validate if endpoint surface has reached
        if (Math.Abs(phreaticLine.Points.Last().X - maxXCoordinateSurface) > GeometryConstants.Tolerance)
        {
            var endPoint = new PlLinePoint(maxXCoordinateSurface, phreaticLine.Points.Last().Z);
            phreaticLine.Points.Add(endPoint);
        }
    }

    private void AddPhreaticLineAlongSurfaceLevel(PlLine phreaticLine)
    {
        // Add phreatic point at 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.
        int surfacePointIndex = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder)) + 1;
        var intersected = false;
        for (int i = surfacePointIndex; i < SurfaceLine.Geometry.Points.Count; i++)
        {
            double z = Math.Max(WaterLevelPolder, SurfaceLine.Geometry.Points[i].Z - PlLineOffsetBelowDikeToeAtPolder);
            double x = SurfaceLine.Geometry.Points[i].X;
            if (WaterLevelPolder > SurfaceLine.Geometry.Points[i].Z - PlLineOffsetBelowDikeToeAtPolder)
            {
                // Determine intersection between would be phreatic segment and polderlevel. Add that as next point.
                var waterLevelPolderLine = new Line<Point2D>(new Point2D(SurfaceLine.Geometry.Points.First().X, WaterLevelPolder),
                                                             new Point2D(SurfaceLine.Geometry.Points.Last().X, WaterLevelPolder));
                var slopeLine = new Line<Point2D>(new Point2D(phreaticLine.Points[phreaticLine.Points.Count - 1].X, phreaticLine.Points[phreaticLine.Points.Count - 1].Z),
                                                  new Point2D(SurfaceLine.Geometry.Points[i].X, SurfaceLine.Geometry.Points[i].Z - PlLineOffsetBelowDikeToeAtPolder));
                Point2D intersectionPoint;
                if (waterLevelPolderLine.IntersectsZ(slopeLine, out intersectionPoint))
                {
                    x = intersectionPoint.X;
                }

                intersected = true;
            }

            var point = new PlLinePoint(x, z);
            phreaticLine.Points.Add(point);
            if (intersected)
            {
                break;
            }
        }
    }

    private void AdjustForDitchAtPolderSide(PlLine phreaticLine, int surfacePointIndex)
    {
        const double maxDouble = 99999.999;
        var phreaticPolderPartialLine = new Line();
        //#bka: hier niet langer ook starten met waterlevel als waterlevel onder bottomditch zit!
        phreaticPolderPartialLine.SetBeginAndEndPoints(new Point2D(phreaticLine.Points[0].X, WaterLevelPolder),
                                                       new Point2D(maxDouble, WaterLevelPolder));
        AddIntersectionDitchDikeSegmentPolderLevelToPhreatic(phreaticLine, surfacePointIndex, phreaticPolderPartialLine);
        AddIntersectionDitchPolderSegmentPolderLevelToPhreatic(phreaticLine, phreaticPolderPartialLine);
    }

    /// <summary>
    /// Intersection between two line segments:
    /// -Ditchpolder Surfaceline segment
    /// -Polder waterlevel
    /// </summary>
    /// <param name="phreaticLine"></param>
    /// <param name="phreaticPolderPartialLine"></param>
    private void AddIntersectionDitchPolderSegmentPolderLevelToPhreatic(PlLine phreaticLine, Line phreaticPolderPartialLine)
    {
        GeometryPoint pointDitchPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide);
        if (pointDitchPolder != null)
        {
            int indexatDitchPolder = SurfaceLine.Geometry.Points.IndexOf(pointDitchPolder);
            var lineDitchPolderSide = new Line();
            if (indexatDitchPolder > 1)
            {
                lineDitchPolderSide.SetBeginAndEndPoints(new Point2D(SurfaceLine.Geometry.Points[indexatDitchPolder - 1].X, SurfaceLine.Geometry.Points[indexatDitchPolder - 1].Z),
                                                         new Point2D(SurfaceLine.Geometry.Points[indexatDitchPolder].X, SurfaceLine.Geometry.Points[indexatDitchPolder].Z));

                var intersectDitchPolderPhreatic = new GeometryPoint();
                if (LineHelper.DetermineStrictIntersectionPoint(lineDitchPolderSide, phreaticPolderPartialLine, ref intersectDitchPolderPhreatic))
                {
                    phreaticLine.Points.Add(new PlLinePoint(intersectDitchPolderPhreatic.X, intersectDitchPolderPhreatic.Z));
                }
                else
                {
                    phreaticLine.Points.Add(new PlLinePoint(pointDitchPolder.X, phreaticPolderPartialLine.EndPoint.Z));
                }
            }
            else
            {
                throw new PlLinesCreatorException("DitchPolderSide should be part of a line segment");
            }
        }
    }

    /// <summary>
    /// Intersection between two line segments:
    /// -DitchDike Surfaceline segment
    /// -Polder waterlevel
    /// </summary>
    /// <param name="phreaticLine"></param>
    /// <param name="surfacePointIndex"></param>
    /// <param name="phreaticPolderPartialLine"></param>
    private void AddIntersectionDitchDikeSegmentPolderLevelToPhreatic(PlLine phreaticLine, int surfacePointIndex, Line phreaticPolderPartialLine)
    {
        if (surfacePointIndex + 1 < SurfaceLine.Geometry.Points.Count)
        {
            var lineDitchDikeSide = new Line();
            lineDitchDikeSide.SetBeginAndEndPoints(new Point2D(SurfaceLine.Geometry.Points[surfacePointIndex].X, SurfaceLine.Geometry.Points[surfacePointIndex].Z),
                                                   new Point2D(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X, SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z));

            var intersectDitchDikePhreatic = new GeometryPoint();
            if (LineHelper.DetermineStrictIntersectionPoint(lineDitchDikeSide, phreaticPolderPartialLine, ref intersectDitchDikePhreatic))
            {
                phreaticLine.Points.Add(new PlLinePoint(intersectDitchDikePhreatic.X, intersectDitchDikePhreatic.Z));
            }
            else
            {
                phreaticLine.Points.Add(new PlLinePoint(SurfaceLine.Geometry.Points[surfacePointIndex].X, phreaticPolderPartialLine.EndPoint.Z));
            }
        }
        else
        {
            throw new PlLinesCreatorException("Could not create DikeSegmentPolderLevel");
        }
    }

    private PlLinePoint DetermineIntersectionBetweenPolderLevelAndDike(double polderLevel)
    {
        var polderlevelLine = new Line();
        double startXCoordinate = SurfaceLine.Geometry.Points.OrderBy(p => p.X).First().X;
        GeometryPoint pointEndOfprofile = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside);
        polderlevelLine.SetBeginAndEndPoints(new Point2D(startXCoordinate, polderLevel),
                                             new Point2D(pointEndOfprofile.X, polderLevel));

        ThrowWhenWaterLevelAboveDike(polderLevel, SurfaceLine);
        // Find intersection between dike top (riverside because that is taken as default DTH) and dike toe (polder side) with polder level.
        // Start at dike toe and work your way back to dike top.
        int startPosition = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder));
        int endPosition = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver));
        for (int surfacePointIndex = startPosition; surfacePointIndex > endPosition; surfacePointIndex--)
        {
            var surfaceLineSegment = new Line();
            surfaceLineSegment.SetBeginAndEndPoints(new Point2D(SurfaceLine.Geometry.Points[surfacePointIndex - 1].X,
                                                                SurfaceLine.Geometry.Points[surfacePointIndex - 1].Z), new Point2D(SurfaceLine.Geometry.Points[surfacePointIndex].X,
                                                                                                                                   SurfaceLine.Geometry.Points[surfacePointIndex].Z));
            var intersectPoint = new GeometryPoint();
            if (LineHelper.DetermineStrictIntersectionPoint(surfaceLineSegment, polderlevelLine, ref intersectPoint))
            {
                return new PlLinePoint(intersectPoint.X, intersectPoint.Z);
            }
        }

        // When polder level is below all points in the given surface line segment, there is no intersection point.
        return null;
    }

    /// <summary>
    /// Detects whether the phreatic line is above the surface line between DikeTopAtRiver and DikeToeAtPolder
    /// 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 ValidatePhreaticBelowDike(PlLine phreaticLine)
    {
        int startIndex = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver));
        int stopIndex = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder));
        var 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 * SurfaceLine.Geometry.Points.Count);
        var iterationCount = 0;
        while (!stopIteration)
        {
            iterationCount++;
            var foundIntersect = false;
            for (var phreaticPointIndex = 0; phreaticPointIndex < phreaticLine.Points.Count - 1; phreaticPointIndex++)
            {
                var phreaticLineSegment = new Line();
                phreaticLineSegment.SetBeginAndEndPoints(
                    new Point2D(phreaticLine.Points[phreaticPointIndex].X, phreaticLine.Points[phreaticPointIndex].Z),
                    new Point2D(phreaticLine.Points[phreaticPointIndex + 1].X, phreaticLine.Points[phreaticPointIndex + 1].Z));
                for (int surfacePointIndex = startIndex; surfacePointIndex < stopIndex; surfacePointIndex++)
                {
                    var surfaceLineSegment = new Line();
                    surfaceLineSegment.SetBeginAndEndPoints(
                        new Point2D(SurfaceLine.Geometry.Points[surfacePointIndex].X, SurfaceLine.Geometry.Points[surfacePointIndex].Z),
                        new Point2D(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X, SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z));
                    var intersectGeoPoint = new GeometryPoint();
                    var intersectPoint = new Point2D();
                    if (LineHelper.DetermineStrictIntersectionPoint(phreaticLineSegment, surfaceLineSegment, ref intersectGeoPoint))
                    {
                        intersectPoint.X = intersectGeoPoint.X;
                        intersectPoint.Z = intersectGeoPoint.Z;
                        // Prevent any adding when intersectPoint is already on Pl
                        if (!intersectPoint.LocationEquals(phreaticLineSegment.BeginPoint) &&
                            !intersectPoint.LocationEquals(phreaticLineSegment.EndPoint))
                        {
                            // 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 > WaterLevelPolder + cOffsetPhreaticLineBelowSurface)
                            {
                                var newPlLinePoint = new PlLinePoint(intersectPoint.X,
                                                                     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(phreaticPointIndex + 1, newPlLinePoint);
                                    if (SurfaceLine.Geometry.Points[surfacePointIndex + 1].X.IsNearEqual(intersectPoint.X))
                                    {
                                        // phreatic point and surfaceline point are postioned above each other, so replace the phreatic point with the intersect point
                                        phreaticLine.Points[phreaticPointIndex + 2] =
                                            new PlLinePoint(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X,
                                                            SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z -
                                                            cOffsetPhreaticLineBelowSurface);
                                    }
                                    else
                                    {
                                        var newNextPlLinePoint =
                                            new PlLinePoint(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X,
                                                            SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z -
                                                            cOffsetPhreaticLineBelowSurface);
                                        if (newNextPlLinePoint.X <= phreaticLine.Points[phreaticPointIndex + 2].X)
                                        {
                                            if (newNextPlLinePoint.X.IsNearEqual(phreaticLine.Points[phreaticPointIndex + 2].X))
                                            {
                                                // If phreatic point already exist on this x-coordinate just replace it with the new point
                                                phreaticLine.Points[phreaticPointIndex + 2] = newNextPlLinePoint;
                                            }
                                            else
                                            {
                                                phreaticLine.Points.Insert(phreaticPointIndex + 2,
                                                                           newNextPlLinePoint);
                                            }
                                        }
                                    }

                                    foundIntersect = true;
                                }

                                break;
                            }
                        }
                    }
                }
            }

            if (!foundIntersect)
            {
                stopIteration = true;
            }

            if (iterationCount > cMaxIterations)
            {
                throw new PlLinesCreatorException("PlLinesCreator.ValidatePhreaticBelowDike() cannot succeed");
            }
        }
    }

    /// <summary>
    /// Create PL1 (is phreatic level)
    /// </summary>
    /// <returns></returns>
    /// 
    private PlLine CreatePlLine1ByExpertKnowledge()
    {
        ThrowWhenWaterLevelPolderAboveDikeTopAtPolder();

        //Create Phreatic line and add polderwater level 
        double startXCoordinate = SurfaceLine.Geometry.Points.OrderBy(p => p.X).First().X;

        var phreaticLine = new PlLine();
        phreaticLine.IsPhreatic = true;
        double? waterLevel = WaterLevelToUse();
        // Waterlevel MUST be above of Dike Toe At River Side, if not, an error is thrown
        SurfaceLine.CheckWaterLevelIsAboveDikeToeRiverSide(waterLevel);
        phreaticLine.Points.Add(new PlLinePoint(startXCoordinate, waterLevel.Value));

        PlLinePoint intersectionPhreaticDike = null;

        //Determine intersection Phreaticlevel - Berm (or dike)
        intersectionPhreaticDike = DetermineIntersectionBetweenWaterLevelAndDike(waterLevel.Value);
        if (intersectionPhreaticDike == null)
        {
            throw new PlLinesCreatorException("DetermineIntersectionBetweenWaterLevelAndDike failed");
        }

        //Add intersectioncoordinate to phreatic line list if not already in list.
        if (!phreaticLine.Points.Contains(intersectionPhreaticDike))
        {
            phreaticLine.Points.Add(intersectionPhreaticDike);
        }

        //Complete phreatic line                                                       
        if (IsUseLowWaterLevel)
        {
            CreatePhreaticLineSegmentsInsideDikeForLowRiverLevel(phreaticLine, waterLevelRiverLow.Value, WaterLevelRiverHigh);
        }
        else
        {
            CreatePhreaticLineSegmentsInsideDikeForHighRiverLevel(phreaticLine);
        }

        CreatePhreaticLineSegmentsInShoulderAndPolder(phreaticLine);

        //Check if phreatic line is above 
        ValidatePhreaticAboveWaterLevelPolder(phreaticLine);
        //Check if phreatic line is below the surface line
        ValidatePhreaticBelowDike(phreaticLine);

        // currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
        CurrentPl1Line = phreaticLine;
        return phreaticLine;
    }

    /// <summary>
    /// Check if phreatic line is above waterlevel
    /// If it is below, correct it
    /// </summary>
    /// <param name="phreaticLine"></param>
    private void ValidatePhreaticAboveWaterLevelPolder(PlLine phreaticLine)
    {
        foreach (PlLinePoint phreaticPoint in phreaticLine.Points)
        {
            if (phreaticPoint.Z < WaterLevelPolder)
            {
                phreaticPoint.Z = WaterLevelPolder;
            }
        }
    }

    private SoilProfile1D GetRelevantSoilProfileForAquiferLayersSearch()
    {
        GeometryPoint relevantPoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
        if (relevantPoint == null)
        {
            throw new PlLinesCreatorException("Characteristic point \"dike toe at polder\" is not defined.");
        }

        return GetSoilProfileBelowPoint(relevantPoint.X);
    }

    /// <summary>
    /// Throws the exception when water level above dike.
    /// </summary>
    /// <param name="waterLevelRiver">The water level river.</param>
    /// <param name="surfaceLine2">The surface line.</param>
    /// <exception cref="PlLinesCreatorException"></exception>
    private void ThrowWhenWaterLevelAboveDike(double waterLevelRiver, SurfaceLine2 surfaceLine2)
    {
        double dikeTopAtRiverHeight = surfaceLine2.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).Z;
        if (waterLevelRiver > dikeTopAtRiverHeight)
        {
            throw new PlLinesCreatorException($"Phreatic water level should have an intersection with the dike at least once (phreatic line higher ({waterLevelRiver:F2} m) than surface line ({dikeTopAtRiverHeight:F2}))");
        }
    }

    /// <summary>
    /// 
    /// </summary>
    private void ThrowWhenWaterLevelPolderAboveDikeTopAtPolder()
    {
        double dikeTopPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).Z;
        if (WaterLevelPolder > dikeTopPolder)
        {
            throw new PlLinesCreatorException($"Waterlevel ({WaterLevelPolder:0.00}) in polder higher than dike top at polder ({dikeTopPolder:0.00}) in surfaceline {SurfaceLine.Name}");
        }
    }

    private void ValidateSoilProfileForPlLinesCreation()
    {
        var validator = new SoilProfileValidator
        {
            SurfaceLine = SurfaceLine,
            SoilProfileType = SoilProfileType,
            SoilProfile1D = SoilProfile,
            SoilProfile2D = SoilProfile2D,
            DikeEmbankmentMaterial = DikeEmbankmentMaterial
        };

        validator.ValidateSoilProfileForPlLinesCreator();
    }

    internal bool IsSurfaceLineIntersectedByAquiferAtPolder(PlLineType plLineType, out GeometryPoint intersectionPoint)
    {
        // If all the layers below the dike toe at polder are aquifer, no intersection point can be found 
        SoilProfile1D relevantSoilProfile = GetRelevantSoilProfileForAquiferLayersSearch();
        if (relevantSoilProfile.Layers.All(layer => layer.IsAquifer))
        {
            intersectionPoint = null;
            return false;
        }
        
        // Soil profile 2D
        double xStartSearch = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
        if (SoilProfileType == SoilProfileType.ProfileType2D)
        {
            SoilProfile2DHelper.LayerType aquiferType = plLineType == PlLineType.Pl3 ? SoilProfile2DHelper.LayerType.BottomAquiferCluster : SoilProfile2DHelper.LayerType.InBetweenAquiferCluster;
            return SoilProfile2DHelper.IsSurfaceLineIntersectedByAquifer(aquiferType, SoilProfile2D, xStartSearch, out intersectionPoint);
        }
        
        // Soil profile 1D
        SoilLayer1D aquifer = GetRelevantAquiferLayer(plLineType, relevantSoilProfile);
        if (aquifer != null)
        {
            int indexStartSearch = SurfaceLine.Geometry.Points.FindIndex(p => p.X.IsNearEqual(xStartSearch, toleranceAlmostEquals));
            double xIntersection = SurfaceLine.Geometry.GetXatZStartingAt(indexStartSearch, aquifer.TopLevel);
            if (!double.IsNaN(xIntersection))
            {
                intersectionPoint = new GeometryPoint
                {
                    X = xIntersection,
                    Z = aquifer.TopLevel
                };
                return true;
            }
        }

        intersectionPoint = null;
        return false;
    }


    /// <summary>
    /// The intersection of PL 1 with the ditch is found starting the search at <see cref="intersectionDitchAquifer"/>.
    /// If PL 1 (polder level) is higher than <see cref="intersectionDitchAquifer"/>, search in the outward direction.
    /// If PL 1 (polder level) is lower than <see cref="intersectionDitchAquifer"/>, search in the inward direction.
    /// Return the first intersection point found as <see cref="intersectionDitchPl1"/>. 
    /// </summary>
    /// <param name="intersectionDitchAquifer">The intersection point of the ditch with the aquifer.</param>
    /// <param name="intersectionDitchPl1">Returns the intersection point of the ditch with PL 1. If no intersection point is
    /// found, returns null.</param>
    /// <returns><c>true</c> if the ditch intersects the polder level; otherwise, <c>false</c>.</returns>
    internal bool IsDitchIntersectedByPl1(GeometryPoint intersectionDitchAquifer, out GeometryPoint intersectionDitchPl1)
    {
        intersectionDitchPl1 = new GeometryPoint();
        bool isPl1AboveAquifer = CurrentPl1Line.Points.Last().Z.IsGreaterThanOrEqualTo(intersectionDitchAquifer.Z, toleranceAlmostEquals);
        GeometryPoint dikeToeAtPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
        double xStart = SurfaceLine.Geometry.Points.First(p => p.X.IsGreaterThanOrEqualTo(intersectionDitchAquifer.X, toleranceAlmostEquals)).X;
        int indexStart = SurfaceLine.Geometry.Points.FindIndex(p => p.X.AlmostEquals(xStart, toleranceAlmostEquals)) - 1;
        int indexDikeToeAtPolder = SurfaceLine.Geometry.Points.FindIndex(p => p.X.AlmostEquals(dikeToeAtPolder.X, toleranceAlmostEquals));
        int indexEnd = isPl1AboveAquifer ? indexDikeToeAtPolder : SurfaceLine.Geometry.Points.Count - 1;

        if ((isPl1AboveAquifer && indexStart < indexEnd) || (!isPl1AboveAquifer && indexStart > indexEnd))
        {
            throw new PlLinesCreatorException(Resources.PlLinesCreator_ErrorDuringIsDitchIntersectedByPl1);
        }

        int index = indexStart;
        while (index != indexEnd)
        {
            GeometryPoint currentPoint = SurfaceLine.Geometry.Points[index];
            GeometryPoint nextPoint = SurfaceLine.Geometry.Points[index + 1];
            var surfaceLineSegment = new Line();
            surfaceLineSegment.SetBeginAndEndPoints(new Point2D(currentPoint.X, currentPoint.Z), new Point2D(nextPoint.X, nextPoint.Z));
            var phreaticPolderPartialLine = new Line();
            PlLinePoint endPl1 = CurrentPl1Line.Points.First(p => p.X.IsGreaterThanOrEqualTo(currentPoint.X, toleranceAlmostEquals));
            PlLinePoint beginPl1 = CurrentPl1Line.Points.Last(p => p.X.IsLessThan(currentPoint.X, toleranceAlmostEquals));
            phreaticPolderPartialLine.SetBeginAndEndPoints(new Point2D(beginPl1.X, beginPl1.Z), new Point2D(endPl1.X, endPl1.Z));
            if (LineHelper.DetermineStrictIntersectionPoint(surfaceLineSegment, phreaticPolderPartialLine, ref intersectionDitchPl1))
            {
                return true;
            }
            index = isPl1AboveAquifer ? index - 1 : index + 1;
        }

        intersectionDitchPl1 = null;
        return false;
    }
}