// 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.KernelWrappers.Common;
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 CharacteristicPointType = Deltares.DamEngine.Data.Geotechnics.CharacteristicPointType;
using SoilProfile = Deltares.DamEngine.Data.Geotechnics.SoilProfile;

namespace Deltares.DamEngine.Calculators.KernelWrappers.MacroStabilityCommon;

public static class PlLinesToWaternetConverter
{
    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)
    {
        // The validation of the soil profile was already performed in the Pl-lines creator. 
        // Here we only have to check that the PL-lines were indeed created.
        ThrowWhenPlLinesIsNull(plLines);
        
        double xLeft = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside).X;
        double xRight = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X;

        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)
    {
        // The validation of the soil profile was already performed in the Pl-lines creator. 
        // Here we only have to check that the PL-lines were indeed created.
        ThrowWhenPlLinesIsNull(plLines);
        
        SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.BottomAquiferCluster, soilProfile, out Point2D[] bottomAquiferCoordinates, out _);
        
        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;
    }
    
    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 (!SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int inBetweenAquiferCount))
        {
            return;
        }

        var headLine = CreateLine<HeadLine>(plLine4, headLine4Name);
        waternet.HeadLineList.Add(headLine);
        for (var i = 0; i < inBetweenAquiferCount; i++)
        {
            SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] inBetweenAquiferUpperCoordinates, out Point2D[] inBetweenAquiferLowerCoordinates, i);
            
            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));
            }
        }
    }

    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();
        AvoidUpwardsVerticalLine(line.CalcPoints);
        line.SyncPoints();

        return line;
    }
   
    /// <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:
            {
                SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.HighestAquiferCluster, soilProfile2D, out Point2D[] highestAquiferCoordinates, out _);
                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:
            {
                SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.LowestLayer, soilProfile2D, out _, out Point2D[] lowestBoundaryCoordinates);
                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 IsBelowSoilProfile(SoilProfile1D soilProfile, PlLine line)
    {
        double bottomSoilProfileLevel = soilProfile.BottomLevel;
        return line.Points.Any(point => point.Z <= bottomSoilProfileLevel);
    }

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

    /// <summary>
    /// This correction is needed for the creation of the STIX file because a reference line containing a vertical part
    /// is always drawn in downwards direction by D-Stability. 
    /// </summary>
    internal static void AvoidUpwardsVerticalLine(IList<Point2D> points)
    {
        for (var i = 0; i < points.Count - 1; i++)
        {
            if (!points[i].X.IsNearEqual(points[i + 1].X, toleranceAlmostEqual) || points[i].Z.IsGreaterThanOrEqualTo(points[i + 1].Z, toleranceAlmostEqual))
            {
                continue;
            }
            
            if (i == points.Count - 2 || (points[i + 1].X + 0.01).IsLessThan(points[i + 2].X, toleranceAlmostEqual))
            {
                points[i + 1].X += 0.01;
            }
            else
            {
                points[i + 1].X = points[i + 2].X - (points[i + 2].X - points[i + 1].X) / 2;
            }
        }
    }
}