// 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.Data;
using System.Linq;
using Deltares.DamEngine.Calculators.Properties;
using Deltares.DamEngine.Data.General;
using Deltares.DamEngine.Data.General.PlLines;
using Deltares.DamEngine.Data.Geometry;
using Deltares.DamEngine.Data.Geotechnics;
using Deltares.DamEngine.Data.Standard;
using Deltares.MacroStability.CSharpWrapper.Input;
using CharacteristicPointType = Deltares.DamEngine.Data.Geotechnics.CharacteristicPointType;
using SoilProfile = Deltares.DamEngine.Data.Geotechnics.SoilProfile;

namespace Deltares.DamEngine.Calculators.KernelWrappers.MacroStabilityCommon;

public static class PlLinesToWaternetConverter
{
    public enum LayerType
    {
        BottomAquifer,
        TopLayerInBetweenAquiferCluster,
        BottomLayerInBetweenAquiferCluster,
        HighestAquifer,
        LowestLayer
    }

    private const double deviationX = 1e-06;
    private const double toleranceAlmostEqual = 1e-09;
    
    private const string waternetLine1Name = "Waternet line phreatic line";
    private const string waternetLine2Name = "Penetration zone lower aquifer";
    private const string waternetLine3Name = "Waternet line lower aquifer";
    private const string waternetLine4TopName = "Waternet line in-between aquifer top";
    private const string waternetLine4BottomName = "Waternet line in-between aquifer bottom";
    private const string headLine1Name = "Phreatic line (PL 1)";
    private const string headLine2Name = "Head line 2 (PL 2)";
    private const string headLine3Name = "Head line 3 (PL 3)";
    private const string headLine4Name = "Head line 4 (PL 4)";

    /// <summary>
    /// Creates a <see cref="Waternet"/> based on a 1D profile by assigning the  <see cref="PlLines"/> to waternet lines.
    /// </summary>
    /// <param name="plLines">The pl lines.</param>
    /// <param name="soilProfile1D">The 1D soil profile.</param>
    /// <param name="surfaceLine">The surface line.</param>
    /// <param name="penetrationLength">Length of the penetration.</param>
    /// <param name="pressureType">The type of distribution of the vertical water pressures.</param>
    /// <returns>A <see cref="Waternet"/>.</returns>
    /// <exception cref="NoNullAllowedException">Thrown when <paramref name="plLines"/>, <paramref name="soilProfile1D"/> or
    /// <paramref name="surfaceLine"/> is <c>null</c>.</exception>
    public static Waternet CreateWaternetBasedOnPlLines(PlLines plLines, SoilProfile1D soilProfile1D, SurfaceLine2 surfaceLine,
                                                        double penetrationLength, IntrusionVerticalWaterPressureType? pressureType)
    {
        ThrowWhenPlLinesIsNull(plLines);
        ThrowWhenSoilProfileIsNull(soilProfile1D);
        ThrowWhenSurfaceLineIsNull(surfaceLine);

        double xLeft = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside).X;
        double xRight = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X;

        if (soilProfile1D.BottomAquiferLayer == null)
        {
            throw new NoNullAllowedException(string.Format(Resources.NoBottomAquiferLayer, soilProfile1D.Name));
        }

        var waternet = new Waternet();

        PlLine plLine = plLines.Lines[PlLineType.Pl1];
        var headLine = CreateLine<HeadLine>(plLine, headLine1Name);
        if (plLine != null && !IsBelowSoilProfile(soilProfile1D, plLine))
        {
            waternet.PhreaticLine = CreateLine<PhreaticLine>(plLine, headLine1Name);
            WaternetLine waternetLine = CreateWaternetLineForPhreaticLine(soilProfile1D, surfaceLine.Geometry, headLine, xLeft, xRight, pressureType);
            waternetLine.Name = waternetLine1Name;
            waternetLine.HeadLine = headLine;
            waternetLine.HeadLine.SyncCalcPoints();
            waternet.WaternetLineList.Add(waternetLine);
        }

        plLine = plLines.Lines[PlLineType.Pl2];
        headLine = CreateLine<HeadLine>(plLine, headLine2Name);
        if (headLine != null && !IsBelowSoilProfile(soilProfile1D, plLine) &&
            pressureType == IntrusionVerticalWaterPressureType.SemiTimeDependent && penetrationLength > 0)
        {
            waternet.HeadLineList.Add(headLine);
            double level = soilProfile1D.BottomAquiferLayer.TopLevel + penetrationLength;
            WaternetLine waternetLine = CreateWaternetLine(level, xLeft, xRight);
            waternetLine.Name = waternetLine2Name;
            waternetLine.HeadLine = headLine;
            waternet.WaternetLineList.Add(waternetLine);
        }

        plLine = plLines.Lines[PlLineType.Pl3];
        headLine = CreateLine<HeadLine>(plLine, headLine3Name);
        if (headLine != null && !IsBelowSoilProfile(soilProfile1D, plLine) &&
            pressureType != IntrusionVerticalWaterPressureType.FullHydroStatic)
        {
            waternet.HeadLineList.Add(headLine);
            double level = soilProfile1D.BottomAquiferLayer.TopLevel;
            WaternetLine waternetLine = CreateWaternetLine(level, xLeft, xRight);
            waternetLine.Name = waternetLine3Name;
            waternetLine.HeadLine = headLine;
            waternet.WaternetLineList.Add(waternetLine);
        }

        CreateWaternetLinesForInBetweenAquifers(waternet, soilProfile1D, xLeft, xRight, plLines.Lines[PlLineType.Pl4], pressureType);

        AdjustWaternetLineOfPhreaticLineWhenCoincidingWithOtherWaternetLines(waternet);

        return waternet;
    }

    /// <summary>
    /// Creates a <see cref="Waternet"/> based on a 2D profile by assigning the  <see cref="PlLines"/> to waternet lines.
    /// </summary>
    /// <param name="plLines">The <see cref="PlLines"/> to convert.</param>
    /// <param name="soilProfile">The <see cref="SoilProfile2D"/> to convert the <see cref="PlLines"/> with.</param>
    /// <param name="penetrationLength">The penetration length.</param>
    /// <param name="pressureType">The type of distribution of the vertical water pressures.</param>
    /// <returns>A <see cref="Waternet"/>.</returns>
    /// <exception cref="NoNullAllowedException">Thrown when <paramref name="plLines"/> or <paramref name="soilProfile"/>
    /// is <c>null</c>.</exception>
    public static Waternet CreateWaternetBasedOnPlLines(PlLines plLines, SoilProfile2D soilProfile, double penetrationLength, IntrusionVerticalWaterPressureType? pressureType)
    {
        ThrowWhenPlLinesIsNull(plLines);
        ThrowWhenSoilProfileIsNull(soilProfile);

        Point2D[] bottomAquiferCoordinates = DetermineLayerBoundaryPoints(LayerType.BottomAquifer, soilProfile).ToArray();
        if (!bottomAquiferCoordinates.Any())
        {
            throw new NoNullAllowedException(string.Format(Resources.NoBottomAquiferLayer, soilProfile.Name));
        }

        var waternet = new Waternet();
        PlLine plLine = plLines.Lines[PlLineType.Pl1];
        if (plLine != null)
        {
            var headLine = CreateLine<HeadLine>(plLine, headLine1Name);
            waternet.PhreaticLine = CreateLine<PhreaticLine>(plLine, headLine1Name);
            WaternetLine waternetLine = CreateWaternetLineForPhreaticLine(soilProfile, plLine, pressureType);
            waternetLine.Name = waternetLine1Name;
            waternetLine.HeadLine = headLine;
            waternet.WaternetLineList.Add(waternetLine);
        }

        plLine = plLines.Lines[PlLineType.Pl2];
        if (plLine != null && pressureType == IntrusionVerticalWaterPressureType.SemiTimeDependent && penetrationLength > 0)
        {
            var headLine = CreateLine<HeadLine>(plLine, headLine2Name);
            waternet.HeadLineList.Add(headLine);
            WaternetLine waternetLine = CreateWaternetLine(bottomAquiferCoordinates, penetrationLength);
            waternetLine.Name = waternetLine2Name;
            waternetLine.HeadLine = headLine;
            waternet.WaternetLineList.Add(waternetLine);
        }

        plLine = plLines.Lines[PlLineType.Pl3];
        if (plLine != null && pressureType != IntrusionVerticalWaterPressureType.FullHydroStatic)
        {
            var headLine = CreateLine<HeadLine>(plLine, headLine3Name);
            waternet.HeadLineList.Add(headLine);

            WaternetLine waternetLine = CreateWaternetLine(bottomAquiferCoordinates);
            waternetLine.HeadLine = headLine;
            waternetLine.Name = waternetLine3Name;
            waternet.WaternetLineList.Add(waternetLine);
        }

        CreateWaternetLinesForInBetweenAquifers(waternet, soilProfile, plLines.Lines[PlLineType.Pl4], pressureType);

        AdjustWaternetLineOfPhreaticLineWhenCoincidingWithOtherWaternetLines(waternet);

        return waternet;
    }

    internal static WaternetLine CreateWaternetLine(double level, double xLeft, double xRight)
    {
        var waternetLine = new WaternetLine();
        waternetLine.Points.Add(new GeometryPoint(xLeft, level));
        waternetLine.Points.Add(new GeometryPoint(xRight, level));
        waternetLine.SyncCalcPoints();
        return waternetLine;
    }

    /// <summary>
    /// Determine all the xCoordinates to make cross-sections. 
    /// </summary>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <returns>All the xCoordinates of the soil profile 2D.</returns>
    private static double[] DetermineAllXCoordinatesOfSoilProfile(SoilProfile2D soilProfile)
    {
        IEnumerable<Point2D> points = soilProfile.Surfaces.SelectMany(surf => surf.GeometrySurface.OuterLoop.CalcPoints);
        double[] xCoordinates = points.Select(point => point.X).OrderBy(x => x).Distinct().ToArray();
        return xCoordinates;
    }

    private static void CreateWaternetLinesForInBetweenAquifers(Waternet waternet, SoilProfile1D soilProfile1D, double xLeft, double xRight, PlLine plLine4, IntrusionVerticalWaterPressureType? pressureType)
    {
        var headLine = CreateLine<HeadLine>(plLine4, headLine4Name);
        if (headLine == null || IsBelowSoilProfile(soilProfile1D, plLine4))
        {
            return;
        }

        if (soilProfile1D.GetInBetweenAquiferClusters == null || pressureType == IntrusionVerticalWaterPressureType.FullHydroStatic)
        {
            return;
        }

        waternet.HeadLineList.Add(headLine);
        var number = 0;
        foreach ((SoilLayer1D, SoilLayer1D) inBetweenAquiferCluster in soilProfile1D.GetInBetweenAquiferClusters.Where(layer => layer is { Item1: not null, Item2: not null }))
        {
            number++;
            double levelTop = inBetweenAquiferCluster.Item1.TopLevel;
            WaternetLine waternetLineTop = CreateWaternetLine(levelTop, xLeft, xRight);
            waternetLineTop.Name = waternetLine4TopName + " (" + number + ")";
            waternetLineTop.HeadLine = headLine;
            waternet.WaternetLineList.Add(waternetLineTop);
            
            double levelBottom = inBetweenAquiferCluster.Item2.BottomLevel;
            WaternetLine waternetLineBottom = CreateWaternetLine(levelBottom, xLeft, xRight);
            waternetLineBottom.Name = waternetLine4BottomName + " (" + number + ")";
            waternetLineBottom.HeadLine = headLine;
            waternet.WaternetLineList.Add(waternetLineBottom);
        }
    }

    private static void CreateWaternetLinesForInBetweenAquifers(Waternet waternet, SoilProfile2D soilProfile, PlLine plLine4,
                                                                IntrusionVerticalWaterPressureType? pressureType)
    {
        if (plLine4 == null)
        {
            return;
        }

        if (pressureType == IntrusionVerticalWaterPressureType.FullHydroStatic)
        {
            return;
        }

        if (!AreInBetweenAquiferClustersPresent(soilProfile, out int inBetweenAquiferCount))
        {
            return;
        }

        var headLine = CreateLine<HeadLine>(plLine4, headLine4Name);
        waternet.HeadLineList.Add(headLine);
        for (var i = 0; i < inBetweenAquiferCount; i++)
        {
            Point2D[] inBetweenAquiferUpperCoordinates = DetermineLayerBoundaryPoints(LayerType.TopLayerInBetweenAquiferCluster, soilProfile, i).ToArray();
            Point2D[] inBetweenAquiferLowerCoordinates = DetermineLayerBoundaryPoints(LayerType.BottomLayerInBetweenAquiferCluster, soilProfile, i).ToArray();

            if (inBetweenAquiferUpperCoordinates.Any() && inBetweenAquiferLowerCoordinates.Any())
            {
                WaternetLine waternetLineTop = CreateWaternetLine(inBetweenAquiferUpperCoordinates);
                waternetLineTop.HeadLine = headLine;
                waternetLineTop.Name = waternetLine4TopName + " (" + (i + 1) + ")";
                waternet.WaternetLineList.Add(waternetLineTop);
                
                WaternetLine waternetLineBottom = CreateWaternetLine(inBetweenAquiferLowerCoordinates);
                waternetLineBottom.HeadLine = headLine;
                waternetLineBottom.Name = waternetLine4BottomName + " (" + (i + 1) + ")";
                waternet.WaternetLineList.Add(waternetLineBottom);
            }
            else
            {
                throw new NoNullAllowedException(string.Format(Resources.PlLinesToWaternetConverter_ErrorDuringDeterminationCoordinatesInBetweenAquifer, soilProfile.Name));
            }
        }
    }

    internal static bool AreInBetweenAquiferClustersPresent(SoilProfile2D soilProfile, out int count)
    {
        double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile);
        var currentCount = 0;
        var previousCount = 0;
        for (var i = 0; i < xCoordinates.Length; i++)
        {
            SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinates[i]);
            currentCount = crossSection.GetInBetweenAquiferClusters?.Count ?? 0;
            if (i > 0 && currentCount != previousCount)
            {
                count = 0;
                return false;
            }

            previousCount = currentCount;
        }

        count = currentCount;
        return count > 0;
    }

    private static TLineType CreateLine<TLineType>(PlLine plLine, string name)
        where TLineType : GeometryPointString, new()
    {
        if (plLine == null || !plLine.Points.Any())
        {
            return null;
        }

        var line = new TLineType();
        line.Points.AddRange(plLine.Points);
        line.SyncCalcPoints();
        line.Name = name;
        return line;
    }

    private static WaternetLine CreateWaternetLine(IEnumerable<Point2D> coordinates, double zOffSet = 0)
    {
        var line = new WaternetLine();
        foreach (Point2D coordinate in coordinates)
        {
            var point = new GeometryPoint(coordinate.X, coordinate.Z + zOffSet);
            line.Points.Add(point);
        }

        line.SyncCalcPoints();
        LineHelper.RemoveDuplicatedPoints(line.CalcPoints, toleranceAlmostEqual);
        line.RemoveUnnecessaryPoints();
        line.SyncPoints();

        return line;
    }

    /// <summary>
    /// Determines all the points of the layer boundary for the specified layer type:
    /// - For LayerType = BottomAquifer, it will return the top boundary of the bottom aquifer (= waternet line for PL3)
    /// - For TopLayerInBetweenAquiferCluster, it will return the top boundary of the in-between aquifer (= waternet line for PL4)
    /// - For BottomLayerInBetweenAquiferCluster, it will return the bottom boundary of the in-between aquifer (= waternet line for PL4)
    /// - For HighestAquifer, it will return the top boundary of the highest aquifer (= waternet line for PL1 for "Hydrostatic" waternet)
    /// - For LowestLayer, it will return the bottom boundary of the lowest layer (= waternet line for PL1 for "Full hydrostatic" waternet)
    /// </summary>
    /// <param name="layerType">The type of layer boundary.</param>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <param name="indexInBetweenAquifer">In case of several in-between aquifers, the index of the in-between aquifer must be specified.</param>
    /// <returns>All the points of the layer boundary.</returns>
    internal static IEnumerable<Point2D> DetermineLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, int indexInBetweenAquifer = -1)
    {
        double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile);
        SoilLayer1D previousAquiferLayer = null;
        var xCoordinatesAll = new List<double>();
        foreach (double xCoordinate in xCoordinates)
        {
            if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile, indexInBetweenAquifer))
            {
                xCoordinatesAll.Add(xCoordinate - deviationX);
                xCoordinatesAll.Add(xCoordinate + deviationX);
            }
            else
            {
                xCoordinatesAll.Add(xCoordinate);
            }
        }

        var layerBoundaryPoints = new List<Point2D>();
        foreach (double xCoordinate in xCoordinatesAll)
        {
            SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate);

            // Determine if the layer is in range of the previous layer
            // If not, return empty coordinates, because the layer is interrupted
            SoilLayer1D currentAquifer = GetSoilLayer1D(layerType, crossSection, indexInBetweenAquifer);
            if (previousAquiferLayer != null && currentAquifer != null
                                             && !AreHorizontallyConnected(previousAquiferLayer, currentAquifer))
            {
                return Enumerable.Empty<Point2D>();
            }

            if (currentAquifer != null)
            {
                previousAquiferLayer = currentAquifer;
                layerBoundaryPoints.Add(layerType is LayerType.BottomLayerInBetweenAquiferCluster or LayerType.LowestLayer ? new Point2D(xCoordinate, currentAquifer.BottomLevel) : new Point2D(xCoordinate, currentAquifer.TopLevel));
            }
        }

        // Perform a short validation that the coordinates are fully defined from the beginning to the end 
        // of the profile. If not, the layer is not continuous.
        if (!IsLayerBoundaryContinuous(layerBoundaryPoints, xCoordinates.First(), xCoordinates.Last()))
        {
            return Enumerable.Empty<Point2D>();
        }

        return layerBoundaryPoints;
    }

    private static bool IsLayerBoundaryContinuous(List<Point2D> boundaryPoints, double leftGeometryBoundary, double rightGeometryBoundary)
    {
        return boundaryPoints.Any() && (Math.Abs(boundaryPoints.First().X - leftGeometryBoundary) < toleranceAlmostEqual)
                                    && Math.Abs(boundaryPoints.Last().X - rightGeometryBoundary) < toleranceAlmostEqual;
    }

    private static SoilLayer1D GetSoilLayer1D(LayerType layerType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer)
    {
        switch (layerType)
        {
            case LayerType.BottomAquifer:
                return soilProfile1D.BottomAquiferLayer;
            case LayerType.TopLayerInBetweenAquiferCluster:
                return soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item1;
            case LayerType.BottomLayerInBetweenAquiferCluster:
                return soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item2;
            case LayerType.HighestAquifer:
                return soilProfile1D.GetHighestAquifer();
            case LayerType.LowestLayer:
                return soilProfile1D.Layers.Last();
        }

        return null;
    }

    /// <summary>
    /// The position of the waternet line associated to the phreatic line depends on the water pressure type:
    /// - For DAM Standard type, this line lies on the bottom of the soil layers “in which the phreatic plane lies”.
    /// - For Linear type, it is the surface line.
    /// - For Hydrostatic type, it is the top of the highest aquifer.
    /// - For Full hydrostatic, it is the bottom of the lowest layer.
    /// - For Semi time dependent, idem DAM Standard.
    /// </summary>
    /// <param name="soilProfile1D"></param>
    /// <param name="surfaceLine"></param>
    /// <param name="phreaticLine"></param>
    /// <param name="xLeft"></param>
    /// <param name="xRight"></param>
    /// <param name="pressureType"></param>
    /// <returns></returns>
    private static WaternetLine CreateWaternetLineForPhreaticLine(SoilProfile1D soilProfile1D, GeometryPointString surfaceLine,
                                                                  HeadLine phreaticLine, double xLeft, double xRight,
                                                                  IntrusionVerticalWaterPressureType? pressureType)
    {
        switch (pressureType)
        {
            case IntrusionVerticalWaterPressureType.Linear:
            {
                List<Point2D> coordinates = surfaceLine.Points.Select(point => new Point2D(point.X, point.Z)).ToList();
                return CreateWaternetLine(coordinates);
            }
            case IntrusionVerticalWaterPressureType.HydroStatic:
            {
                double level = soilProfile1D.GetHighestAquifer().TopLevel;
                return CreateWaternetLine(level, xLeft, xRight);
            }
            case IntrusionVerticalWaterPressureType.FullHydroStatic:
            {
                double level = soilProfile1D.Layers.Last().BottomLevel;
                return CreateWaternetLine(level, xLeft, xRight);
            }
            case IntrusionVerticalWaterPressureType.Standard:
            case IntrusionVerticalWaterPressureType.SemiTimeDependent:
            {
                double level = soilProfile1D.GetLayerAt(phreaticLine.GetMinZ()).BottomLevel;
                // Make sure the waternet line for PL 1 is not below the one for PL 3
                double bottomAquiferLevel = soilProfile1D.BottomAquiferLayer.TopLevel;
                level = Math.Max(level, bottomAquiferLevel);
                return CreateWaternetLine(level, xLeft, xRight);
            }
            case null:
                break;
            default:
                throw new ArgumentOutOfRangeException(nameof(pressureType), pressureType, null);
        }

        return null;
    }

    /// <summary>
    /// The position of the waternet line associated to the phreatic line depends on the water pressure type:
    /// - For DAM Standard type, this line lies on the bottom of the soil layers “in which the phreatic plane lies”.
    /// - For Linear type, it is the surface line.
    /// - For Hydrostatic type, it is the top of the highest aquifer.
    /// - For Full hydrostatic, it is the bottom of the lowest layer.
    /// - For Semi time dependent, idem DAM Standard.
    /// </summary>
    /// <param name="soilProfile2D"></param>
    /// <param name="plLine"></param>
    /// <param name="pressureType"></param>
    /// <returns></returns>
    /// <exception cref="ArgumentOutOfRangeException"></exception>
    private static WaternetLine CreateWaternetLineForPhreaticLine(SoilProfile2D soilProfile2D, PlLine plLine,
                                                                  IntrusionVerticalWaterPressureType? pressureType)
    {
        switch (pressureType)
        {
            case IntrusionVerticalWaterPressureType.Standard:
            case IntrusionVerticalWaterPressureType.SemiTimeDependent:
            {
                return CreateWaternetLineForPhreaticLineForDamStandard(soilProfile2D, plLine);
            }
            case IntrusionVerticalWaterPressureType.HydroStatic:
            {
                Point2D[] highestAquiferCoordinates = DetermineLayerBoundaryPoints(LayerType.HighestAquifer, soilProfile2D).ToArray();
                return CreateWaternetLine(highestAquiferCoordinates);
            }
            case IntrusionVerticalWaterPressureType.Linear:
            {
                List<Point2D> xCoordinates = soilProfile2D.Geometry.SurfaceLine.Points.Select(point => new Point2D(point.X, point.Z)).ToList();
                return CreateWaternetLine(xCoordinates);
            }
            case IntrusionVerticalWaterPressureType.FullHydroStatic:
            {
                Point2D[] lowestBoundaryCoordinates = DetermineLayerBoundaryPoints(LayerType.LowestLayer, soilProfile2D).ToArray();
                return CreateWaternetLine(lowestBoundaryCoordinates);
            }
            case null:
                break;
            default:
                throw new ArgumentOutOfRangeException(nameof(pressureType), pressureType, null);
        }

        return null;
    }

    /// <summary>
    /// Create the waternet line associated to the phreatic line (PL1) for the waternet type "Dam Standard".
    /// This line lies on the bottom of the soil layers “in which the phreatic plane lies”.
    /// Where PL1 lies above the ground level, the line lies on the ground level.
    /// As a consequence, the line can contain vertical "jumps" when layers have vertical parts.
    /// </summary>
    /// <param name="soilProfile2D"></param>
    /// <param name="plLine"></param>
    /// <returns></returns>
    internal static WaternetLine CreateWaternetLineForPhreaticLineForDamStandard(SoilProfile2D soilProfile2D, PlLine plLine)
    {
        var waternetLine = new WaternetLine();
        var curves = new List<GeometryPointString>();
        foreach (SoilLayer2D surface in soilProfile2D.Surfaces)
        {
            GeometryLoop outerLoop = surface.GeometrySurface.OuterLoop;
            // Add the first point to get a closed GeometryPointString
            outerLoop.CalcPoints.Add(new Point2D(outerLoop.CalcPoints[0].X, outerLoop.CalcPoints[0].Z));

            var plLineInGeometryPoint = new GeometryPointString();
            foreach (PlLinePoint point in plLine.Points)
            {
                plLineInGeometryPoint.Points.Add(new GeometryPoint(point.X, point.Z));
                plLineInGeometryPoint.SyncCalcPoints();
            }

            List<Point2D> intersectionPoints = outerLoop.IntersectionXzPointsWithGeometryString(plLineInGeometryPoint);
            // Only the layers that intersect the phreatic line are kept.
            if (intersectionPoints.Count > 0)
            {
                var geometrySurface = new GeometrySurface(outerLoop);
                GeometryPointString bottom = geometrySurface.DetermineBottomGeometrySurface();
                bottom.SyncPoints();
                bottom.SortPointsByXAscending();
                curves.Add(bottom);
            }
            
            // Add also the surface line in case of free water
            curves.Add(soilProfile2D.Geometry.SurfaceLine);
        }

        // The waternet line is the polyline formed with the lowest lines connected to each other.
        List<Line> lines = DivideCurvesIntoLines(curves);
        IEnumerable<Point2D> allPoints = curves.SelectMany(curve => curve.CalcPoints);
        double[] xCoordinates = allPoints.Select(point => point.X).OrderBy(x => x).Distinct().ToArray();
        List<Line> splitLines = SplitLinesAtXCoordinates(xCoordinates, lines);
        for (var i = 0; i < xCoordinates.Length - 1; i++)
        {
            List<Line> relevantLinesAtX = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < toleranceAlmostEqual && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0);
            if (relevantLinesAtX.Count == 0)
            {
                continue;
            }
            double minZBegin = relevantLinesAtX.Select(line => line.BeginPoint.Z).Min();
            List<Line> relevantLinesAtXWithZBeginMax = relevantLinesAtX.FindAll(line => Math.Abs(line.BeginPoint.Z - minZBegin) < toleranceAlmostEqual);
            double minZEnd = relevantLinesAtXWithZBeginMax.Select(line => line.EndPoint.Z).Min();
            Line relevantLine = relevantLinesAtXWithZBeginMax.Find(line => Math.Abs(line.EndPoint.Z - minZEnd) < toleranceAlmostEqual);
            waternetLine.Points.Add(new GeometryPoint(relevantLine.BeginPoint.X, relevantLine.BeginPoint.Z));
            waternetLine.Points.Add(new GeometryPoint(relevantLine.EndPoint.X, relevantLine.EndPoint.Z));
        }

        // Make sure the waternet line for PL 1 is not in the bottom aquifer.
        foreach (GeometryPoint waternetLinePoint in waternetLine.Points)
        {
            double bottomAquiferLevel = soilProfile2D.GetSoilProfile1D(waternetLinePoint.X).BottomAquiferLayer.TopLevel;
            waternetLinePoint.Z = Math.Max(waternetLinePoint.Z, bottomAquiferLevel);
        }

        waternetLine.SyncCalcPoints();
        LineHelper.RemoveDuplicatedPoints(waternetLine.CalcPoints, toleranceAlmostEqual);
        waternetLine.RemoveUnnecessaryPoints();
        waternetLine.SyncPoints();

        return waternetLine;
    }

    private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile, int indexInBetweenAquifer)
    {
        SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX);
        SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX);
        if (crossSectionLeft != null && crossSectionRight != null)
        {
            SoilLayer1D currentAquiferLeft = GetSoilLayer1D(layerType, crossSectionLeft, indexInBetweenAquifer);
            SoilLayer1D currentAquiferRight = GetSoilLayer1D(layerType, crossSectionRight, indexInBetweenAquifer);
            if (currentAquiferLeft != null && currentAquiferRight != null)
            {
                if (layerType is LayerType.BottomAquifer or LayerType.TopLayerInBetweenAquiferCluster)
                {
                    return Math.Abs(currentAquiferLeft.TopLevel - currentAquiferRight.TopLevel) > GeometryConstants.Accuracy;
                }

                return Math.Abs(currentAquiferLeft.BottomLevel - currentAquiferRight.BottomLevel) > GeometryConstants.Accuracy;
            }
        }

        return false;
    }

    private static bool AreHorizontallyConnected(SoilLayer1D leftSoilLayer, SoilLayer1D rightSoilLayer)
    {
        // The layers are identical
        if (ReferenceEquals(leftSoilLayer, rightSoilLayer))
        {
            return true;
        }

        // Left soil layer envelopes whole right soil layer
        if (leftSoilLayer.BottomLevel <= rightSoilLayer.BottomLevel
            && leftSoilLayer.TopLevel >= rightSoilLayer.TopLevel)
        {
            return true;
        }

        // Right soil layer envelopes whole left soil layer
        if (rightSoilLayer.BottomLevel <= leftSoilLayer.BottomLevel
            && rightSoilLayer.TopLevel >= leftSoilLayer.TopLevel)
        {
            return true;
        }

        return (rightSoilLayer.TopLevel <= leftSoilLayer.TopLevel && rightSoilLayer.TopLevel >= leftSoilLayer.BottomLevel) // Top level lies inbetween the left soil layer
               || (rightSoilLayer.BottomLevel >= leftSoilLayer.BottomLevel && rightSoilLayer.BottomLevel <= leftSoilLayer.TopLevel); // Bottom level lies inbetween the left soil layer
    }

    private static bool IsBelowSoilProfile(SoilProfile1D soilProfile, PlLine line)
    {
        double bottomSoilProfileLevel = soilProfile.BottomLevel;
        return line.Points.Any(point => point.Z <= bottomSoilProfileLevel);
    }

    /// <summary>
    /// Throws when the soil profile is not assigned.
    /// </summary>
    /// <param name="soilProfile">The soil profile.</param>
    /// <exception cref="NoNullAllowedException"></exception>
    private static void ThrowWhenSoilProfileIsNull(SoilProfile soilProfile)
    {
        if (soilProfile == null)
        {
            throw new NoNullAllowedException(Resources.PlLinesToWaternetConverter_NoSoilProfileDefined);
        }
    }

    /// <summary>
    /// Throws when the pl lines object is not assigned.
    /// </summary>
    /// <param name="plLines">The pl lines.</param>
    /// <exception cref="NoNullAllowedException"></exception>
    private static void ThrowWhenPlLinesIsNull(PlLines plLines)
    {
        if (plLines == null)
        {
            throw new NoNullAllowedException(Resources.PlLinesToWaternetConverter_NoPlLinesDefined);
        }
    }

    /// <summary>
    /// Throws when the surface line object is not assigned.
    /// </summary>
    /// <param name="surfaceLine">The surface line.</param>
    /// <exception cref="NoNullAllowedException"></exception>
    private static void ThrowWhenSurfaceLineIsNull(SurfaceLine2 surfaceLine)
    {
        if (surfaceLine == null)
        {
            throw new NoNullAllowedException(Resources.PlLinesToWaternetConverter_NoSurfaceLineDefined);
        }
    }

    /// <summary>
    /// Waternet lines can't coincide otherwise the program does not know which headline to use.
    /// That's why when the waternet line of the phreatic line coincides with the waternet line of either PL2, PL3 or PL4,
    /// an adjustment of the waternet line of the phreatic line is performed by moving it 1 mm upwards. 
    /// </summary>
    /// <param name="waternet">The waternet.</param>
    private static void AdjustWaternetLineOfPhreaticLineWhenCoincidingWithOtherWaternetLines(Waternet waternet)
    {
        const double minimumDistance = 0.001;
        WaternetLine waternetLine1 = waternet.WaternetLineList.Find(w => w.Name == waternetLine1Name);
        if (waternetLine1 != null)
        {
            foreach (GeometryPoint waternetLine1Point
                     in from waternetLineOther in waternet.WaternetLineList.FindAll(w => w.Name != waternetLine1Name)
                        from waternetLine1Point in waternetLine1.Points
                        where waternetLineOther.Points.Any(point => Math.Abs(point.Z - waternetLine1Point.Z) < toleranceAlmostEqual)
                        select waternetLine1Point)
            {
                waternetLine1Point.Z += minimumDistance;
            }
            waternetLine1.SyncCalcPoints();
        }
    }

    private static List<Line> SplitLinesAtXCoordinates(double[] xCoordinates, List<Line> lines)
    {
        var splitLines = new List<Line>();
        foreach (Line line in lines)
        {
            List<double> xCoordinatesToBeAdded = xCoordinates.Where(xCoordinate => line.BeginPoint.X.IsLessThan(xCoordinate, toleranceAlmostEqual) && xCoordinate.IsLessThan(line.EndPoint.X, toleranceAlmostEqual)).ToList();
            
            if (xCoordinatesToBeAdded.Count > 0)
            {
                SplitLineAtXCoordinate(xCoordinatesToBeAdded, line, splitLines);
            }
            else
            {
                splitLines.Add(line);
            }
        }

        return splitLines;
    }

    internal static void SplitLineAtXCoordinate(List<double> xCoordinatesToBeAdded, Line line, List<Line> splitLines)
    {
        for (var i = 0; i < xCoordinatesToBeAdded.Count; i++)
        {
            var point1 = new Point2D(xCoordinatesToBeAdded[i], Math.Max(line.BeginPoint.Z, line.EndPoint.Z) + 1);
            var point2 = new Point2D(xCoordinatesToBeAdded[i], Math.Min(line.BeginPoint.Z, line.EndPoint.Z) - 1);
            double zCoordinate = line.GetIntersectPointXz(new Line(point1, point2)).Z;

            var middlePoint = new Point2D(xCoordinatesToBeAdded[i], zCoordinate);
            Point2D beginPoint = i == 0 ? line.BeginPoint : splitLines.Last().EndPoint;

            splitLines.Add(new Line(beginPoint, middlePoint));
            if (i == xCoordinatesToBeAdded.Count - 1)
            {
                splitLines.Add(new Line(middlePoint, line.EndPoint));
            }
        }
    }

    private static List<Line> DivideCurvesIntoLines(List<GeometryPointString> curves)
    {
        var lines = new List<Line>();
        foreach (IList<GeometryPoint> points in curves.Select(curve => curve.Points))
        {
            for (var i = 0; i < points.Count - 1; i++)
            {
                var line = new Line
                {
                    BeginPoint = new Point2D(points[i].X, points[i].Z),
                    EndPoint = new Point2D(points[i + 1].X, points[i + 1].Z)
                };
                lines.Add(line);
            }
        }

        return lines;
    }
}