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

using System;
using System.Collections.Generic;
using System.Linq;
using Deltares.Geometry;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Mathematics;
using Deltares.Standard.EventPublisher;
using Deltares.Standard.Extensions;
using Deltares.Standard.Language;
using Deltares.Standard.Logging;
using Deltares.Standard.Validation;

namespace Deltares.Geotechnics.WaternetCreator
{
    /// <summary>
    /// Automatic creation of the waternet according to method Ringtoets WTI 2017.
    /// Refer to the memo of Alexander van Duinen (1209434-012-GEO-0002-m-Schematisering waterspanningen in WTI 2017 (Ringtoets).docx) for more information:
    /// https://repos.deltares.nl/repos/FailureMechanisms/FailureMechanisms/DikesMacroStability/trunk/doc/Work
    /// </summary>
    public class WaternetCreatorWTI : IWaternetCreator
    {
        private Soil bottomAquifer;
        private SoilLayer1D bottomAquiferLayer;
        private Soil inBetweenAquifer;
        private SoilLayer1D inBetweenAquiferLayerBottom;
        private SoilLayer1D inBetweenAquiferLayerTop;
        private Location location;
        private readonly PlLinesCreatorWTI plLinesCreator = new PlLinesCreatorWTI();
        private SoilProfile1D profile1D;
        private SoilProfile2D profile2D;
        private Waternet waternet;
        bool isPl4Created = false;

        /// <summary>
        /// Create a waternet based on a soil profile 2D, a surface line and the output from PL-lines creator, for different scenario's:
        /// Sand dike on sand has only phreatic line.
        /// Clay dike on sand have only a PL-line 2 on the top of first aquifer.
        /// The 2 other scenario's may have an in-between aquifer (PL-line 4) and an in-between aquitard (PL-line 2)
        /// </summary>
        /// <param name="waternet">The actual waternet</param>
        /// <param name="location">The input parameters for waternet creator (Location class)</param>
        /// <returns>The updated waternet</returns>
        public void UpdateWaternet(Waternet waternet, Location location)
        {
            if (!HasWaterLevel(location))
            {
                LogManager.Add(new LogMessage(LogMessageType.Warning, this, "Cannot update waternet if no water level is available."));
                return;
            }

            this.waternet = waternet;
            Dictionary<WaternetLine, string> assignedHeadLines = new Dictionary<WaternetLine, string>();

            // Send events, used for undo manager
            if (waternet.PhreaticLine != null)
            {
                DataEventPublisher.BeforeChange(waternet.PhreaticLine.Points);
                waternet.PhreaticLine.Points.Clear();
            }

            DataEventPublisher.BeforeChange(this.waternet.WaternetLineList);
            DataEventPublisher.BeforeChange(this.waternet.HeadLineList);

            this.waternet = waternet;
            this.waternet.IsGenerated = location.WaternetCreationMode == WaternetCreationMode.CreateWaternet;

            if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
            {
                this.waternet.WaternetLineList.Clear();
            }
            else
            {
                foreach (WaternetLine waternetLine in this.waternet.WaternetLineList)
                {
                    if (waternetLine.HeadLine != null)
                    {
                        assignedHeadLines[waternetLine] = waternetLine.HeadLine.Name;
                    }
                }
            }

            this.waternet.HeadLineList.Clear();
            plLinesCreator.Inwards = location.Inwards;
            plLinesCreator.UnitWeightWater = waternet.UnitWeight;
            AssignLocationParameters(location);
            profile1D = null;
            double xUsedToCreate1DProfile;

            if (location.SoilProfile2D != null)
            {
                profile2D = location.SoilProfile2D;
                GeometryPoint dikeTopRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
                GeometryPoint dikeTopPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);

                if (dikeTopRiver != null && dikeTopPolder != null)
                {
                    xUsedToCreate1DProfile = (dikeTopRiver.X + dikeTopPolder.X) / 2;
                }
                else
                {
                    xUsedToCreate1DProfile = location.Surfaceline.Geometry.GetMinX();
                }
                profile1D = location.SoilProfile2D.GetSoilProfile1D(xUsedToCreate1DProfile);
            }
            else
            {
                profile1D = location.SoilProfile1D;
            }

            // Create PL1 = Phreatic Line and add the surface as a waternet line 
            CreatePhreaticLineAndAssignItToSurfaceLine();

            // For "Sand dike on sand" (2B) only the phreatic line is present, so the waternet creation is finished
            if (location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand)
            {
                if (DetermineAquitards().Any())
                {
                    LogManager.Add(new LogMessage(LogMessageType.Warning, location, LocalizationManager.GetTranslatedText(typeof(WaternetCreator), "SandDikeOnSandWarning")));
                }
            }
            else
            {
                // If no pleistoceen is found, the waternet creation is finished
                DetermineAquifers(profile1D);

                if (bottomAquiferLayer != null)
                {
                    // For scenarios's 1A, 2A and 1B, Head 3 is created and assigned to the Pleistoceen.
                    // However, if no aquitards are present, Head 3 is not created
                    var aquitard = DetermineAquitards();
                    if (aquitard.Any())
                    {
                        // For "Clay dike on sand (1B)", only 1 penetration zone (related to Head 2) is created at the bottom of the toppest aquitard. 
                        // Nothing happens with the other aquitards below. For "Clay/Sand dike on clay (1A and 2A), 2 extra penetration zones (related to Head 2)
                        // are created if an in-between aquifer is present.
                        CreatePl2AndAssignItToPenetrationZones();

                        if (!double.IsNaN(location.HeadInPLLine3))
                        {
                            CreatePl3AndAssignItToDeepestAquifer();

                            // Only for "Clay/Sand dike on clay (1A and 2A), the in-between aquifer (if present) is related to Head 4
                            if (inBetweenAquiferLayerTop != null && inBetweenAquiferLayerBottom != null && !double.IsNaN(location.HeadInPLLine4))
                            {
                                if (location.DikeSoilScenario != DikeSoilScenario.ClayDikeOnSand)
                                {
                                    CreatePl4AndAssignItToInBetweenAquifer();
                                    isPl4Created = true;
                                }
                            }
                        }
                    }
                }
                EnsureHydrostaticPorePressureInLayersWhereThePhreaticLinePasses();
            }

            // Remove head lines which are not used
            if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
            {
                foreach (HeadLine headLine in waternet.HeadLineList.ToArray())
                {
                    if (waternet.WaternetLineList.All(p => p.HeadLine == null || p.HeadLine.Name != headLine.Name))
                    {
                        waternet.HeadLineList.Remove(headLine);
                    }
                }
            }

            // Restore assigned head lines
            if (location.WaternetCreationMode == WaternetCreationMode.CreatePhreaticAndHeadLinesOnly)
            {
                foreach (WaternetLine waternetLine in assignedHeadLines.Keys)
                {
                    if (!waternet.HeadLineList.Contains(waternetLine.HeadLine) && waternetLine.HeadLine != waternet.PhreaticLine)
                    {
                        waternetLine.HeadLine = waternet.HeadLineList.FirstOrDefault(p => p.Name == assignedHeadLines[waternetLine]);
                    }
                }
            }

            DataEventPublisher.AfterChange(waternet, "PhreaticLine");
            DataEventPublisher.DataListModified(waternet.PhreaticLine.Points);

            DataEventPublisher.DataListModified(this.waternet.WaternetLineList);
            DataEventPublisher.DataListModified(this.waternet.HeadLineList);
        }

        /// <summary>
        /// Determines whether the waternet can be generated at the specified location <paramref name="location"/>.
        /// </summary>
        /// <param name="location">The location.</param>
        /// <returns>Returns <c>true</c> if the waternet can be generated at the specified location <paramref name="location"/>; otherwise <c>flase</c>.</returns>
        public bool CanGenerateWaternet(Location location)
        {
            return location.Validate().Where(v => v.MessageType == ValidationResultType.Error).ToArray().Length == 0;
        }

        private bool HasWaterLevel(Location location)
        {
            return !double.IsNaN(location.WaterLevelRiver);
        }

        private void AssignLocationParameters(Location location)
        {
            this.location = location;
            this.location.Surfaceline.SortPoints();
        }

        /// <summary>
        /// Creates the phreatic line (PL 1) and assigns it to the surface line.
        /// </summary>
        private void CreatePhreaticLineAndAssignItToSurfaceLine()
        {
            const string phreaticLineName = "Phreatic Line";
            const string surfaceLineName = "Surface line";

            GeometryPointString str = plLinesCreator.CreatePLLine(PLLineType.PL1, location);

            if (waternet.PhreaticLine == null)
            {
                waternet.PhreaticLine = new PhreaticLine { Name = phreaticLineName };
            }

            waternet.PhreaticLine.Points.AddRange(str.Points);

            WaternetLine line = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == surfaceLineName);
            if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
            {
                line = new WaternetLine { Name = surfaceLineName };

                // Clone points for surface line, because they don't have the same behaviour (enabled)
                for (int i = 0; i < location.Surfaceline.Geometry.Points.Count; i++)
                {
                    var point = location.Surfaceline.Geometry.Points[i];
                    if (!double.IsNaN(point.X))
                    {
                        line.Points.Add((GeometryPoint)point.Clone());
                    }
                }

                waternet.AddWaternetLine(line);
            }

            if (line != null)
            {
                line.HeadLine = waternet.PhreaticLine;
            }
        }

        /// <summary>
        /// Ensures an hydrostatic pore pressure distribution in the layers where the phreatic line passes. 
        /// To do this, extra waternet lines are eventually created by:
        /// - First, creating a list of all layers where the phretic line passes;
        /// - Then, creating a waternet line at the bottom of each listed layer only if no waternet line already exists.
        /// </summary>
        private void EnsureHydrostaticPorePressureInLayersWhereThePhreaticLinePasses()
        {
            var listOfLayersWherePhreaticLinePasses = GetListOfLayersWherePhreaticLinePasses();

            GeometryPointString alreadyExistingWaternetLine;
            if (isPl4Created)
            {
                alreadyExistingWaternetLine = GetContourOfInBetweenClusterOfAquifers();
            }
            else
            {
                alreadyExistingWaternetLine = GetAquiferGeometrySurface(true, GeometryTopBottom.Top);
            }

            List<GeometryPointString> listOfCreatedLines = new List<GeometryPointString>();
            if (listOfLayersWherePhreaticLinePasses.Count > 0)
            {
                foreach (var layer in listOfLayersWherePhreaticLinePasses)
                {
                    GeometryPointString bottomSurfaceOfCurrentLayer = layer.GeometrySurface.DetermineBottomGeometrySurface();
                    bool coincide = BottomOfLayerCoincidesWithTheRelevantExistingWaternetLine(layer, alreadyExistingWaternetLine);
                    if (!coincide)
                    {
                        listOfCreatedLines.Add(bottomSurfaceOfCurrentLayer);
                    }
                }
            }

            // if no waternet line already exists, a waternet line is created for each line listed in listOfCreatedLines. 
            if (listOfCreatedLines.Count != 0)
            {
                int number = 1;

                foreach (var line in listOfCreatedLines)
                {
                    string bottomHydrostaticName = "Bottom hydrostatic zone";

                    if (listOfCreatedLines.Count > 1)
                    {
                        bottomHydrostaticName = "Bottom hydrostatic zone (" + number + ")";
                    }
                    number++;

                    WaternetLine waternetLine = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == bottomHydrostaticName);
                    if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                    {
                        waternetLine = new WaternetLine { Name = bottomHydrostaticName };

                        // Clone points 
                        for (int i = 0; i < line.Points.Count; i++)
                        {
                            var point = line.Points[i];
                            if (!double.IsNaN(point.X))
                            {
                                waternetLine.Points.Add((GeometryPoint)point.Clone());
                            }
                        }

                        waternet.AddWaternetLine(waternetLine);
                    }

                    if (waternetLine != null)
                    {
                        waternetLine.HeadLine = waternet.PhreaticLine;
                    }
                }
            }
        }

        private List<SoilLayer2D> GetListOfLayersWherePhreaticLinePasses()
        {
            List<SoilLayer2D> listOfLayersWherePhreaticLinePasses = new List<SoilLayer2D>();

            // First 

            if (waternet.PhreaticLine.Points.Count >= 2)
            {
                for (int i = 0; i < waternet.PhreaticLine.Points.Count - 1; i++)
                {
                    Point2D phreaticLinePoint1 = new Point2D
                    {
                        X = waternet.PhreaticLine.Points[i].X,
                        Y = waternet.PhreaticLine.Points[i].Z
                    };
                    Point2D phreaticLinePoint2 = new Point2D
                    {
                        X = waternet.PhreaticLine.Points[i + 1].X,
                        Y = waternet.PhreaticLine.Points[i + 1].Z
                    };

                    foreach (var layer in profile2D.Surfaces)
                    {
                        for (int j = 0; j < layer.GeometrySurface.OuterLoop.Points.Count - 1; j++)
                        {
                            Point2D layerPoint1 = new Point2D
                            {
                                X = layer.GeometrySurface.OuterLoop.Points[j].X,
                                Y = layer.GeometrySurface.OuterLoop.Points[j].Z
                            };
                            Point2D layerPoint2 = new Point2D
                            {
                                X = layer.GeometrySurface.OuterLoop.Points[j + 1].X,
                                Y = layer.GeometrySurface.OuterLoop.Points[j + 1].Z
                            };
                            var ip = new Point2D();
                            var res = Routines2D.DetermineIf2DLinesIntersectStrickly(phreaticLinePoint1, phreaticLinePoint2, layerPoint1, layerPoint2, ref ip);
                            bool linesIntersect = (res == LineIntersection.Intersects);
                            if (linesIntersect && !listOfLayersWherePhreaticLinePasses.Contains(layer))
                            {
                                listOfLayersWherePhreaticLinePasses.Add(layer);
                            }
                        }
                    }
                }
            }
            return listOfLayersWherePhreaticLinePasses;
        }

        /// <summary>
        /// Determines wheter the bottom of layer <see cref="layer"/> coincides with the relevant existing  waternet line <see cref="alreadyExistingWaternetLine"/>.
        /// </summary>
        /// <param name="layer">The layer.</param>
        /// <param name="alreadyExistingWaternetLine">The already existing waternet line.</param>
        /// <returns><c>true</c> if the bottom of <see cref="layer"/> coincides with the relevant existing waternet line
        /// <see cref="alreadyExistingWaternetLine"/>. Otherwise returns <c>false</c>.</returns>
        private bool BottomOfLayerCoincidesWithTheRelevantExistingWaternetLine(SoilLayer2D layer, GeometryPointString alreadyExistingWaternetLine)
        {
            GeometryPointString bottomSurfaceOfCurrentLayer = layer.GeometrySurface.DetermineBottomGeometrySurface();
            var numberOfCommonPoints = 0;

            // Check if surface bottomSurfaceOfCurrentLayer as common points with the existingWaternelLine. 
            if (alreadyExistingWaternetLine != null)
            {
                for (int j = 0; j < bottomSurfaceOfCurrentLayer.Count; j++)
                {
                    var pointBottomSurface = bottomSurfaceOfCurrentLayer[j];
                    for (int k = 0; k < alreadyExistingWaternetLine.Count; k++)
                    {
                        var pointWaternetLine = alreadyExistingWaternetLine[k];
                        if ((Math.Abs(pointWaternetLine.X - pointBottomSurface.X) < GeometryConstants.Accuracy)
                            && (Math.Abs(pointWaternetLine.Y - pointBottomSurface.Y) < GeometryConstants.Accuracy)
                            && (Math.Abs(pointWaternetLine.Z - pointBottomSurface.Z) < GeometryConstants.Accuracy))
                        {
                            numberOfCommonPoints++;
                        }
                    }
                }
            }
            if (numberOfCommonPoints >= 2)
            {
                return true;
            }
            else
            {
                return false;
            }

        }

        /// <summary>
        /// Create PL2 (i.e. Head 2) and:
        /// - For "Clay dike on clay (1A)" and "Sand dike on clay (2B)", create a penetration zone (with nr. 1) above the Pleistocene
        /// and assign headline Head 2 to this penetration zone. If an in-between cluster of aquifers is present, 2 extra penetration zones are created 
        /// (nr. 2 and 3) around the in-between cluster of aquifers.
        /// - For "Clay dike on sand (1B)", create only 1 penetration zone (with nr. 1) above the toppest aquifer layer (aquifer layers situated 
        /// in the dike are not relevant).
        /// </summary> 
        private void CreatePl2AndAssignItToPenetrationZones()
        {
            const string penetrationZone1Name = "Penetration zone 1";
            const string penetrationZone2Name = "Penetration zone 2";
            const string penetrationZone3Name = "Penetration zone 3";

            double penetrationLength = location.PenetrationLength;
            GeometryPointString pl2 = plLinesCreator.CreatePLLine(PLLineType.PL2, location);

            if (pl2.Points.Count > 0)
            {
                GeometryPointString aquiferGeometryTopBelowPenetrationZone1 = GetAquiferGeometrySurface(true, GeometryTopBottom.Top);

                WaternetLine line2 = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == penetrationZone1Name);
                if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                {
                    line2 = new WaternetLine { Name = penetrationZone1Name };
                    if (aquiferGeometryTopBelowPenetrationZone1 != null)
                    {
                        line2.Points.AddRange(aquiferGeometryTopBelowPenetrationZone1.Points);

                        foreach (var geometryPoint in line2.Points)
                        {
                            geometryPoint.Z = geometryPoint.Z + penetrationLength;
                        }
                        RemovePlLinePointsAboveSurfaceLine(line2);
                    }

                    waternet.AddWaternetLine(line2);
                }

                HeadLine headLine2 = null;

                if (aquiferGeometryTopBelowPenetrationZone1 != null)
                {
                    headLine2 = new HeadLine { Name = "Head 2" };
                    headLine2.Points.Clear();
                    headLine2.Points.AddRange(pl2.Points);

                    waternet.HeadLineList.Add(headLine2);

                    if (line2 != null)
                    {
                        line2.HeadLine = headLine2;
                    }
                }

                // For "Clay dike on clay (1A)" and "Sand dike on clay (2B)", if an in-between cluster of aquifers is present, 
                // 2 extra penetration zones are created (nr. 2 and 3) around the in-between cluster of aquifers.
                if (location.DikeSoilScenario != DikeSoilScenario.ClayDikeOnSand)
                {
                    if (inBetweenAquiferLayerTop != null && inBetweenAquiferLayerBottom != null)
                    {
                        var inBetweenAquiferGeometryTop = GetAquiferGeometrySurface(false, GeometryTopBottom.Top);
                        var inBetweenAquiferGeometryBottom = GetAquiferGeometrySurface(false, GeometryTopBottom.Bottom);

                        WaternetLine line2B = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == penetrationZone3Name);
                        if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                        {
                            line2B = new WaternetLine { Name = penetrationZone3Name };

                            line2B.Points.AddRange(inBetweenAquiferGeometryTop.Points);
                            foreach (var geometryPoint in line2B.Points)
                            {
                                geometryPoint.Z = geometryPoint.Z + penetrationLength;
                            }
                            RemovePlLinePointsAboveSurfaceLine(line2B);

                            waternet.AddWaternetLine(line2B);
                        }

                        if (line2B != null)
                        {
                            line2B.HeadLine = headLine2;
                        }

                        WaternetLine line2C = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == penetrationZone2Name);
                        if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                        {
                            line2C = new WaternetLine { Name = penetrationZone2Name };

                            line2C.Points.AddRange(inBetweenAquiferGeometryBottom.Points);
                            foreach (var geometryPoint in line2C.Points)
                            {
                                geometryPoint.Z = geometryPoint.Z - penetrationLength;
                            }

                            waternet.AddWaternetLine(line2C);
                        }

                        if (line2C != null)
                        {
                            line2C.HeadLine = headLine2;
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Remove the points situated "in the air" (i.e. above surface line) for the current waternet line 
        /// </summary>
        /// <param name="waternetLine">the current waternet line</param>
        private void RemovePlLinePointsAboveSurfaceLine(WaternetLine waternetLine)
        {
            var intersectionPointsBetweenPenetrationZone1AndSurfaceLine =
                DetermineIntersectionPointsBetweenWaternetLineAndSurfaceLine(location.Surfaceline, waternetLine);
            if (intersectionPointsBetweenPenetrationZone1AndSurfaceLine != null)
            {
                foreach (var geometryPoint in intersectionPointsBetweenPenetrationZone1AndSurfaceLine)
                {
                    waternetLine.Points.Add(geometryPoint);
                }
                waternetLine.Points.Sort();
                for (int i = 0; i < waternetLine.Count; i++)
                {
                    var zSurfaceLine = location.Surfaceline.Geometry.GetZAtX(waternetLine.Points[i].X);
                    // if Z co-ordinate are very close, this means that it is an intersection point, so the point is not removed
                    if ((Math.Abs(waternetLine.Points[i].Z - zSurfaceLine) > 0.00000001) && (waternetLine.Points[i].Z > zSurfaceLine))
                    {
                        waternetLine.Points.Remove(waternetLine.Points[i]);
                        i--;
                    }
                }
            }
        }

        /// <summary>
        /// Create PL3 (i.e. Head 3) and assign it to the top surface of the deepest cluster of aquifers 
        /// </summary>
        private void CreatePl3AndAssignItToDeepestAquifer()
        {
            const string pl3Name = "pl3";

            GeometryPointString pl3 = plLinesCreator.CreatePLLine(PLLineType.PL3, location);

            if (pl3 != null)
            {
                WaternetLine line3 = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == pl3Name);
                if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                {
                    line3 = new WaternetLine { Name = pl3Name };

                    // The waternet corresonding to Head 3 is always the top surface of the deepest cluster of aquifers
                    GeometryPointString bottomAquiferGeometryTop = GetAquiferGeometrySurface(true, GeometryTopBottom.Top);

                    line3.Points.AddRange(bottomAquiferGeometryTop.Points);

                    waternet.AddWaternetLine(line3);
                }

                HeadLine headLine3 = new HeadLine { Name = "Head 3" };
                headLine3.Points.Clear();
                headLine3.Points.AddRange(pl3.Points);

                waternet.HeadLineList.Add(headLine3);

                if (line3 != null)
                {
                    line3.HeadLine = headLine3;
                }
            }
        }

        /// <summary>
        /// Create PL4 (i.e. Head 4) and assign it to the 'in between aquifer' when present.
        /// </summary>
        private void CreatePl4AndAssignItToInBetweenAquifer()
        {
            const string pl4Name = "In between aquifer";
            const string pl4BName = "In between aquifer B";

            GeometryPointString pl4 = plLinesCreator.CreatePLLine(PLLineType.PL4, location);
            if (pl4 != null)
            {
                if (pl4.Points.Count > 0)
                {
                    WaternetLine line4 = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == pl4Name);
                    if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                    {
                        line4 = new WaternetLine { Name = pl4Name };

                        if (profile2D != null)
                        {
                            var points = new List<GeometryPoint>();
                            var contour = GetContourOfInBetweenClusterOfAquifers();
                            points.AddRange(contour.Points);
                            line4.Points.AddRange(points);
                        }
                        else
                        {
                            var zLevelTop = GetZLevelAquifer1DProfile(GeometryTopBottom.Top, false);
                            var pointstring = CreateGeometryPointString(zLevelTop);
                            line4.Points.AddRange(pointstring.Points);
                        }

                        waternet.AddWaternetLine(line4);
                    }

                    HeadLine headLine4 = new HeadLine { Name = "Head 4" };
                    headLine4.Points.Clear();
                    headLine4.Points.AddRange(pl4.Points);
                    waternet.HeadLineList.Add(headLine4);

                    if (line4 != null)
                    {
                        line4.HeadLine = headLine4;
                    }

                    if (profile2D == null)
                    {
                        WaternetLine line4B = this.waternet.WaternetLineList.FirstOrDefault(p => p.Name == pl4BName);
                        if (location.WaternetCreationMode == WaternetCreationMode.CreateWaternet)
                        {
                            line4B = new WaternetLine { Name = pl4BName };

                            var zLevelBottom = GetZLevelAquifer1DProfile(GeometryTopBottom.Bottom, false);
                            var pointstringB = CreateGeometryPointString(zLevelBottom);
                            line4B.Points.AddRange(pointstringB.Points);

                            waternet.AddWaternetLine(line4B);
                        }

                        if (line4B != null)
                        {
                            line4B.HeadLine = headLine4;
                        }
                    }
                }
            }
        }

        /// <summary>
        /// If IsBottomAquifer = true, get the top or bottom geometry surface of the toppest aquifer in the deep cluster of aquifers
        /// If IsBottomAquifer = false, get the top or bottom geometry surface of the toppest aquifer in the in-between cluster of aquifers
        /// </summary>
        /// <param name="isBottomAquifer">Choose between true (the toppest aquifer in the deep cluster of aquifers) or false 
        /// (toppest aquifer in the in-between cluster of aquifers)</param>
        /// <param name="topBottom">Choose between GeometryTopBottom.Top and GeometryTopBottom.Bottom</param>
        /// <returns>The geometry surface (top or bottom) of the selected aquifer (bottom or in-between)</returns>
        private GeometryPointString GetAquiferGeometrySurface(bool isBottomAquifer, GeometryTopBottom topBottom)
        {
            if (profile2D != null)
            {
                GeometrySurface aquiferSurface;
                if (isBottomAquifer)
                {
                    aquiferSurface = GeometrySurfaceAquifer(true, true);
                }
                else
                {
                    if (topBottom == GeometryTopBottom.Top)
                    {
                        aquiferSurface = GeometrySurfaceAquifer(false, true);
                    }
                    else
                    {
                        aquiferSurface = GeometrySurfaceAquifer(false, false);
                    }
                }
                if (aquiferSurface != null)
                {
                    GeometryPointString result;
                    if (topBottom == GeometryTopBottom.Top)
                    {
                        result = aquiferSurface.DetermineTopGeometrySurface();
                    }
                    else
                    {
                        result = aquiferSurface.DetermineBottomGeometrySurface();
                    }
                    return result;
                }
                else
                {
                    return null;
                }
            }
            else
            {
                var zLevel = GetZLevelAquifer1DProfile(topBottom, isBottomAquifer);
                return CreateGeometryPointString(zLevel);
            }
        }

        /// <summary>
        /// Get the contour of the highest cluster of in-between aquifers, used to create the waternet line PL4 associated to Head 4.
        /// Note that an aquifer is always a complete layer (form left to right geometry boundary), 
        /// not an isolated layer (see ValidateAquifers in Location.cs).
        /// </summary>
        /// <returns>List of geometry points (closed polygone)</returns>
        private GeometryPointString GetContourOfInBetweenClusterOfAquifers()
        {
            GeometryPointString result;
            if (profile2D != null)
            {
                GeometrySurface geometrySurfaceTopLayer = GeometrySurfaceAquifer(false, true);
                GeometrySurface geometrySurfaceBottomLayer = GeometrySurfaceAquifer(false, false);

                if (inBetweenAquiferLayerTop != null && inBetweenAquiferLayerBottom != null)
                {
                    GeometryPointString pointsTop = geometrySurfaceTopLayer.DetermineTopGeometrySurface();
                    GeometryPointString pointsBottom = geometrySurfaceBottomLayer.DetermineBottomGeometrySurface();
                    pointsTop.SortPointsByXAscending();
                    pointsBottom.Points.Sort();
                    result = pointsTop;
                    pointsBottom.Points.Sort();
                    for (int i = pointsBottom.Count - 1; i >= 0; i--)
                    {
                        result.Points.Add(pointsBottom[i]);
                    }
                    // add again the first point of the top surface to close the contour
                    result.Points.Add(pointsTop[0]);
                    return result;
                }
                else
                {
                    return null;
                }
            }
            else
            {
                var zLevelTop = GetZLevelAquifer1DProfile(GeometryTopBottom.Top, false);
                var zLevelBottom = GetZLevelAquifer1DProfile(GeometryTopBottom.Bottom, false);
                GeometryPointString pointsTop = CreateGeometryPointString(zLevelTop);
                GeometryPointString pointsBottom = CreateGeometryPointString(zLevelBottom);
                result = pointsTop;
                foreach (var geometryPoint in pointsBottom.Points)
                {
                    result.Points.Add(geometryPoint);
                }
                return result;
            }
        }

        /// <summary>
        /// Get the top surface of the highest cluster of in-between aquifers.
        /// </summary>
        /// <returns>List of geometry points (surface)</returns>
        private GeometryPointString GetTopSurfaceOfInBetweenClusterOfAquifers()
        {
            GeometryPointString result;
            if (profile2D != null)
            {
                GeometrySurface geometrySurfaceTopLayer = GeometrySurfaceAquifer(false, true);
                if (inBetweenAquiferLayerTop != null && inBetweenAquiferLayerBottom != null)
                {
                    GeometryPointString pointsTop = geometrySurfaceTopLayer.DetermineTopGeometrySurface();
                    pointsTop.SortPointsByXAscending();
                    return pointsTop;
                }
                else
                {
                    return null;
                }
            }
            else
            {
                var zLevelTop = GetZLevelAquifer1DProfile(GeometryTopBottom.Top, false);
                GeometryPointString pointsTop = CreateGeometryPointString(zLevelTop);
                result = pointsTop;
                return result;
            }
        }

        /// <summary>
        /// If isBottomAquifer = true, get the top or bottom level of the toppest aquifer in the deep cluster of aquifers
        /// If isBottomAquifer = false, get the top or bottom level of the toppest aquifer in the in-between cluster of aquifers
        /// </summary>
        /// <param name="topBottom">Choose between GeometryTopBottom.Top and GeometryTopBottom.Bottom</param>
        /// <param name="isBottomAquifer">Choose between true (the toppest aquifer in the deep cluster of aquifers) or false 
        /// (toppest aquifer in the in-between cluster of aquifers)</param>
        /// <returns>The level of the selected aquifer (bottom or in-between)</returns>
        private double GetZLevelAquifer1DProfile(GeometryTopBottom topBottom, bool isBottomAquifer)
        {
            double zLevel;
            if (isBottomAquifer)
            {
                if (topBottom == GeometryTopBottom.Top)
                {
                    zLevel = profile1D.Layers.Where(l => l.Soil == bottomAquifer)
                                      .Select(l => l.TopLevel)
                                      .FirstOrDefault();
                }
                else
                {
                    zLevel = profile1D.Layers.Where(l => l.Soil == bottomAquifer)
                                      .Select(l => l.BottomLevel)
                                      .FirstOrDefault();
                }
            }
            else
            {
                if (topBottom == GeometryTopBottom.Top)
                {
                    zLevel =
                        profile1D.Layers.Where(l => l.Soil == inBetweenAquifer)
                                 .Select(l => l.TopLevel)
                                 .FirstOrDefault();
                }
                else
                {
                    zLevel =
                        profile1D.Layers.Where(l => l.Soil == inBetweenAquifer)
                                 .Select(l => l.BottomLevel)
                                 .FirstOrDefault();
                }
            }
            return zLevel;
        }

        private GeometryPointString CreateGeometryPointString(double zLevel)
        {
            var pointString = new GeometryPointString();
            var maxX = location.Surfaceline.Geometry.Points.Max(p => p.X);
            var minX = location.Surfaceline.Geometry.Points.Min(p => p.X);
            pointString.Points.Add(new GeometryPoint(minX, 0, zLevel));
            pointString.Points.Add(new GeometryPoint(maxX, 0, zLevel));
            return pointString;
        }

        /// <summary>
        /// Get the geometry surface of:
        /// - If <paramref name="isBottomAquifer"/> = true, the highest aquifer in the deep cluster of aquifers
        /// - If <paramref name="isBottomAquifer"/> = false, the toppest aquifer in the in-between cluster of aquifers
        /// using the 1D profile generated at the middle of Dike Top.
        /// </summary>
        /// <param name="isBottomAquifer">Choose between true (the toppest aquifer in the deep cluster of aquifers) or false
        /// (toppest aquifer in the in-between cluster of aquifers)</param>
        /// <param name="isTopInBetweenAquifer">if set to <c>true</c> [is top in between aquifer].</param>
        /// <returns>
        /// The geometry surface of the selected aquifer (bottom or in-between)
        /// </returns>
        /// <exception cref="System.Exception">No soil profile is present.</exception>
        private GeometrySurface GeometrySurfaceAquifer(bool isBottomAquifer, bool isTopInBetweenAquifer)
        {
            SoilLayer1D relevantAquiferLayerIn1DProfile = null;
            double xDikeTopAtPolder = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).X;
            double xDikeTopAtRiver = location.Surfaceline.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X;
            double xUsed = (xDikeTopAtRiver + xDikeTopAtPolder) / 2;
            SoilProfile1D soilProfile1DBelowDikeTop = location.SoilProfile2D != null ? location.SoilProfile2D.GetSoilProfile1D(xUsed) : location.SoilProfile1D;

            if (soilProfile1DBelowDikeTop == null)
            {
                throw new System.Exception("No soil profile is present.");
            }
            if (isBottomAquifer)
            {
                // Deepest cluster of sand (Pleistocene)
                relevantAquiferLayerIn1DProfile = soilProfile1DBelowDikeTop.BottomAquiferLayer;
            }
            else
            {
                // Intermediate cluster of sand layers
                if (isTopInBetweenAquifer)
                {
                    relevantAquiferLayerIn1DProfile = soilProfile1DBelowDikeTop.InBetweenAquiferLayer;
                }
                else
                {
                    relevantAquiferLayerIn1DProfile = soilProfile1DBelowDikeTop.BottomLayerOfInBetweenAquiferCluster;
                }
            }

            if (relevantAquiferLayerIn1DProfile != null)
            {
                var middlePoint = new GeometryPoint();
                middlePoint.X = xDikeTopAtPolder;
                middlePoint.Z = (relevantAquiferLayerIn1DProfile.TopLevel + relevantAquiferLayerIn1DProfile.BottomLevel)/2;
                var relevantAquiferLayerIn2DProfile = profile2D.GetSoilLayer(middlePoint);

                return relevantAquiferLayerIn2DProfile.GeometrySurface;
            }
            else
            {
                return null;
            }
        }

        /// <summary>
        /// Determine present aquifers within soil profile 1D / 2D.
        /// If several in-between aquifer clusters are present, the highest cluster is used: InBetweenAquifer is the highest aquifer of the selected cluster and
        /// BottomLayerOfInBetweenAquiferCluster is the lowest aquifer of the selected cluster.
        /// If the deep cluster of aquifers contains several aquifers, the highest aquifer is used as BottomAquifer.
        /// </summary>
        /// <param name="currentProfile1D">The soil profile 1D used to determine the InBetweenAquifer and the BottomAquifer</param>
        private void DetermineAquifers(SoilProfile1D currentProfile1D)
        {
            bottomAquiferLayer = currentProfile1D.BottomAquiferLayer;
            if (bottomAquiferLayer == null)
            {
                return;
            }
            inBetweenAquiferLayerTop = currentProfile1D.InBetweenAquiferLayer;
            inBetweenAquiferLayerBottom = currentProfile1D.BottomLayerOfInBetweenAquiferCluster;
        }

        /// <summary>
        /// determine present aquitards within soil profile 1D / 2D
        /// </summary>
        private List<Soil> DetermineAquitards()
        {
            var aquitards = new List<Soil>();
            if (profile2D != null)
            {
                //IsAquifer is defined on SoilLayer now, not on Soil
                aquitards.AddRange(profile2D.Surfaces.Where(s => s.IsAquifer == false).Select(s => s.Soil));

                // aquitards = profile2D.Surfaces.Where(s => s.IsAquifer == false).Select(s => s.Soil);
            }
            else
            {
                aquitards.AddRange(profile1D.Layers.Where(s => s.IsAquifer == false).Select(s => s.Soil));
            }
            return aquitards;
        }

        /// <summary>
        /// Finds the intersection point of a waternet line and a surface line.
        /// </summary>
        /// <param name="surfaceLine">The evaluated surfaceline.</param>
        /// <param name="waternetLine">The waternet line.</param>
        /// <returns>The list of the intersection points, or null in case no intersection was found.</returns>
        private static List<GeometryPoint> DetermineIntersectionPointsBetweenWaternetLineAndSurfaceLine(SurfaceLine2 surfaceLine, WaternetLine waternetLine)
        {
            var listSurfaceLinePoints = surfaceLine.Geometry.Points;
            var listOfIntersectionPointsFound = new List<GeometryPoint>();
            for (int i = 1; i < listSurfaceLinePoints.Count; i++)
            {
                var surfaceLineSegment = new Line();
                surfaceLineSegment.SetBeginAndEndPoints(
                    new GeometryPoint(listSurfaceLinePoints[i - 1].X, 0, listSurfaceLinePoints[i - 1].Z),
                    new GeometryPoint(listSurfaceLinePoints[i].X, 0, listSurfaceLinePoints[i].Z));
                for (int j = 1; j < waternetLine.Count; j++)
                {
                    var waternetLineSegment = new Line();
                    waternetLineSegment.SetBeginAndEndPoints(
                        new GeometryPoint(waternetLine[j - 1].X, 0, waternetLine[j - 1].Z),
                        new GeometryPoint(waternetLine[j].X, 0, waternetLine[j].Z));
                    var intersectPoint = surfaceLineSegment.GetIntersectPointXZ(waternetLineSegment);
                    if (intersectPoint != null)
                    {
                        listOfIntersectionPointsFound.Add(intersectPoint);
                    }
                }
            }
            if (listOfIntersectionPointsFound.Count != 0)
            {
                return listOfIntersectionPointsFound;
            }
            return null;
        }
    }
}