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

namespace Deltares.Geotechnics.WaternetCreator
{
    public enum GeometryTopBottom
    {
        Top,
        Bottom
    }

    public class WaternetCreator : IWaternetCreator
    {
        private Soil bottomAquifer;
        private SoilLayer1D bottomAquiferLayer;
        private SoilLayer1D inBetweenAquiferLayerTop;
        private SoilLayer1D inBetweenAquiferLayerBottom;
        private GeometryPointString contourOfInBetweenCluster;
        private double geometryXStart;
        private Soil inBetweenAquifer;
        private Location location;
        private PlLinesCreator plLinesCreator = new PlLinesCreator();
        private SoilProfile1D profile1D;
        private SoilProfile2D profile2D;
        private Dictionary<double, double> upliftOutputValuesPl3;
        private Dictionary<double, double> upliftOutputValuesPl4;
        private Waternet waternet;

        public WaternetCreator() { }

        /// <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>
        /// <param name="inwards">Is the inwards orientation</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;

            // 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 = true;
            this.waternet.WaternetLineList.Clear();
            this.waternet.HeadLineList.Clear();
            plLinesCreator.Inwards = location.Inwards;
            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 
            AssignPl1AsPhreaticLine();

            // For "Sand dike on sand" (2B) only the phreatic line is present, so the waternet creation is finished
            if (location.DikeSoilScenario == DikeSoilScenario.SandDikeOnSand)
            {
                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, PL3 is not created
                    var aquitard = DetermineAquitards();
                    if (!double.IsNaN(location.HeadInPLLine3) && aquitard.Any())
                    {
                        AssignPl3();
                    }

                    // To continue, at least one aquitard is needed
                    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.
                        if (aquitard.Any())
                        {
                            AssignPl2();
                        }

                        // 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)
                            {
                                AssignPl4();
                            }
                        }
                        if (location.AdjustPl3And4ForUplift)
                        {
                            upliftOutputValuesPl3 = plLinesCreator.GetUpliftOutputValuesPl3();
                            if (location.DikeSoilScenario != DikeSoilScenario.ClayDikeOnSand)
                            {
                                upliftOutputValuesPl4 = plLinesCreator.GetUpliftOutputValuesPl4();
                            }
                        }
                    }
                }
            }

            // Remove head lines which are not used

            foreach (HeadLine headLine in waternet.HeadLineList.ToArray())
            {
                if (waternet.WaternetLineList.All(p => p.HeadLine != headLine))
                {
                    waternet.HeadLineList.Remove(headLine);
                }
            }

            // Send events

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

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

        /// <summary>
        /// order curves clockwise, curve.EndPoint == nextCurve.HeadPoint
        /// </summary>
        /// <param name="firstCurve"></param>
        /// <param name="boundaryCurves"></param>
        /// <returns></returns>
        public List<GeometryCurve> OrderCurves(GeometryCurve firstCurve, List<GeometryCurve> boundaryCurves)
        {
            // TODO: Code clone of GeometrySurface.OrderCurves; Refactoring recommended
            var ordered = new List<GeometryCurve>();
            ordered.Add(firstCurve);
            boundaryCurves.Remove(firstCurve);
            GeometryCurve current = firstCurve;
            do
            {
                var next = GetNextCurve(current, boundaryCurves);
                if (next != null)
                {
                    ordered.Add(next);
                    current = next;
                }
                else
                {
                    next = GetNextCurveReversed(current, boundaryCurves);
                    ordered.Add(next);
                    current = next;
                }
            } while (boundaryCurves.Count > 0);
            return ordered;
        }

        /// <summary>
        /// get next curve so that curve.EndPoint == nextCurve.HeadPoint
        /// </summary>
        /// <param name="curve"></param>
        /// <param name="boundaryCurves"></param>
        /// <returns></returns>
        public GeometryCurve GetNextCurve(GeometryCurve curve, List<GeometryCurve> boundaryCurves)
        {
            // TODO: Code clone of GeometrySurface.GetNextCurve; Refactoring recommended
            foreach (var geometryCurve in boundaryCurves)
            {
                if (curve.EndPoint == geometryCurve.HeadPoint)
                {
                    boundaryCurves.Remove(geometryCurve);
                    return geometryCurve;
                }
            }
            return null;
        }

        /// <summary>
        /// get uplift factors for x, regarding pl4
        /// </summary>
        /// <returns></returns>
        public Dictionary<double, double> GetUpliftOutputValuesPl4()
        {
            return upliftOutputValuesPl4;
        }

        /// <summary>
        /// get uplift factors for x, regarding pl3
        /// </summary>
        /// <returns></returns>
        public Dictionary<double, double> GetUpliftOutputValuesPl3()
        {
            return upliftOutputValuesPl3;
        }

        public bool CanGenerateWaternet(Location location)
        {
            return location.Validate().Where(v => v.MessageType == ValidationResultType.Error).ToArray().Length == 0;
        }

        /// <summary>
        /// create parallel curve ad distance to each original curve 
        /// </summary>
        /// <param name="curve"></param>
        /// <param name="distance"></param>
        /// <returns></returns>
        internal GeometryCurve CreateParallelCurve(GeometryCurve curve, double distance)
        {
            double lX = 0.0;
            double lZ = 0.0;
            double xFraction = 0.0;
            double zFraction = 0.0;
            double lLength = 0.0;
            double lFraction = 0.0;
            double lDistance = 0.0;
            var point1 = curve.HeadPoint;
            var point2 = curve.EndPoint;
            var point3 = new GeometryPoint();
            var point4 = new GeometryPoint();
            if ((point1.X == point2.X) && (point1.Z == point2.Z))
            {
                return new GeometryCurve(point3, point4);
            }
            else if (Math.Abs(distance) > lDistance)
            {
                lX = point2.X - point1.X;
                lZ = point2.Z - point1.Z;
                lLength = Math.Sqrt((lX * lX) + (lZ * lZ));
                lFraction = distance / lLength;
                xFraction = lFraction * lX;
                zFraction = lFraction * lZ;
                point3.X = point1.X - zFraction;
                point4.X = point2.X - zFraction;
                point3.Z = point1.Z + xFraction;
                point4.Z = point2.Z + xFraction;
            }
            else
            {
                point3.X = point1.X;
                point4.X = point2.X;
                point3.Z = point1.Z;
                point4.Z = point2.Z;
            }
            return new GeometryCurve(point3, point4);
        }

        private bool HasWaterLevel(Location location)
        {
            return !double.IsNaN(location.Inwards || location.PlLineCreationMethod == PlLineCreationMethod.RingtoetsWti2017 ? location.WaterLevelRiver : location.WaterLevelRiverLow);
        }

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

        /// <summary>
        /// Create the phreatic line and assign it to the surface line
        /// </summary>
        private void AssignPl1AsPhreaticLine()
        {
            GeometryPointString str = plLinesCreator.CreatePLLine1ByExpertKnowledge(location);

            if (waternet.PhreaticLine == null)
            {
                waternet.PhreaticLine = new PhreaticLine()
                    {
                        Name = "Phreatic Line"
                    };
            }

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

            var line = new WaternetLine()
                {
                    Name = "Surface line"
                };

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

            line.HeadLine = waternet.PhreaticLine;
            waternet.AddWaternetLine(line);
        }

        /// <summary>
        /// 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 AssignPl2()
        {
            double penetrationLength = location.PenetrationLength;
            GeometryPointString pl2 = plLinesCreator.CreatePLLine2ByExpertKnowledge(location);
            List<GeometryPoint> intersectionPointsBetweenPenetrationZone1AndSurfaceLine;

            if (pl2.Points.Count > 0)
            {
                var line2 = new WaternetLine()
                    {
                        Name = "Penetration zone 1"
                    };
                GeometryPointString aquiferGeometryTopBelowPenetrationZone1;
//                if (location.DikeSoilScenario == DikeSoilScenario.ClayDikeOnSand)
//                {
//                    aquiferGeometryTopBelowPenetrationZone1 = GetAquiferGeometrySurface(false, GeometryTopBottom.Top);
//                }
//                else
//                {
                    aquiferGeometryTopBelowPenetrationZone1 = GetAquiferGeometrySurface(true, GeometryTopBottom.Top);
//                }

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

                    foreach (var geometryPoint in line2.Points)
                    {
                        geometryPoint.Z = geometryPoint.Z + penetrationLength;
                    }
                    RemovePlLinePointsAboveSurfaceLine(line2);
                    waternet.AddWaternetLine(line2);
                    waternet.HeadLineList.Add((HeadLine)line2.HeadLine);
                }

                // 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 line2B = new WaternetLine()
                            {
                                Name = "Penetration zone 3",
                                HeadLine = line2.HeadLine
                            };
                        var inBetweenAquiferGeometryTop = GetAquiferGeometrySurface(false, GeometryTopBottom.Top);
                        line2B.Points.AddRange(inBetweenAquiferGeometryTop.Points);
                        foreach (var geometryPoint in line2B.Points)
                        {
                            geometryPoint.Z = geometryPoint.Z + penetrationLength;
                        }
                        RemovePlLinePointsAboveSurfaceLine(line2B);
                        waternet.AddWaternetLine(line2B);

                        var line2C = new WaternetLine()
                            {
                                Name = "Penetration zone 2",
                                HeadLine = line2.HeadLine
                            };
                        var inBetweenAquiferGeometryBottom = GetAquiferGeometrySurface(false, GeometryTopBottom.Bottom);
                        line2C.Points.AddRange(inBetweenAquiferGeometryBottom.Points);
                        foreach (var geometryPoint in line2C.Points)
                        {
                            geometryPoint.Z = geometryPoint.Z - penetrationLength;
                        }
                        waternet.AddWaternetLine(line2C);
                    }
                }
            }
        }

        /// <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>
        /// Assign head line 3 (i.e. PL3) to the top surface of the deepest cluster of aquifers 
        /// </summary>
        private void AssignPl3()
        {
            GeometryPointString pl3 = plLinesCreator.CreatePLLine3ByExpertKnowledge(location);

            if (pl3 != null)
            {
                GeometryPointString bottomAquiferGeometryTop;
                var line3 = new WaternetLine()
                    {
                        Name = "pl3"
                    };
                line3.HeadLine = new HeadLine()
                    {
                        Name = "Head 3"
                    };
                line3.HeadLine.Points.Clear();
                line3.HeadLine.Points.AddRange(pl3.Points);

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

                line3.Points.AddRange(bottomAquiferGeometryTop.Points);

                waternet.AddWaternetLine(line3);
                waternet.HeadLineList.Add((HeadLine)line3.HeadLine);
            }
        }

        /// <summary>
        /// Assign pl4 line to waternet, when in between aquifer is present
        /// </summary>
        private void AssignPl4()
        {
            GeometryPointString pl4 = plLinesCreator.CreatePLLine4ByExpertKnowledge(location);
            if (pl4 != null)
            {
                if (pl4.Points.Count > 0)
                {
                    var line4 = new WaternetLine()
                        {
                            Name = "In between aquifer"
                        };
                    line4.HeadLine = new HeadLine()
                        {
                            Name = "Head 4"
                        };
                    line4.HeadLine.Points.Clear();
                    line4.HeadLine.Points.AddRange(pl4.Points);
                    var points = new List<GeometryPoint>();
                    if (profile2D != null)
                    {
                        var contour = GetContourOfInBetweenClusterOfAquifers();
                        points.AddRange(contour.Points);
                        line4.Points.AddRange(points);
                        waternet.AddWaternetLine(line4);
                        waternet.HeadLineList.Add((HeadLine)line4.HeadLine);
                    }
                    else
                    {
                        var zLevelTop = GetZLevelAquifer1DProfile(GeometryTopBottom.Top, false);
                        var pointstring = CreateGeometryPointString(zLevelTop);
                        line4.Points.AddRange(pointstring.Points);
                        var line4B = new WaternetLine
                            {
                                Name = "In between aquifer",
                                HeadLine = line4.HeadLine
                            };
                        var zLevelBottom = GetZLevelAquifer1DProfile(GeometryTopBottom.Bottom, false);
                        var pointstringB = CreateGeometryPointString(zLevelBottom);
                        line4B.Points.AddRange(pointstringB.Points);
                        waternet.AddWaternetLine(line4);
                        waternet.AddWaternetLine(line4B);
                        waternet.HeadLineList.Add((HeadLine)line4.HeadLine);
                    }
                }
            }
        }

        /// <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)
            {
                GeometryPointString result;
                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)
                {
                    if (topBottom == GeometryTopBottom.Top)
                    {
                        result = aquiferSurface.DetermineTopGeometrySurface();
                    }
                    else
                    {
                        result = aquiferSurface.DetermineBottomGeometrySurface();
                    }
                    return result;
                }
                else
                {
                    return null;
                }
            }
            else
            {
                double zLevel;
                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
        /// </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
            {
                double zLevelTop, zLevelBottom;
                zLevelTop = GetZLevelAquifer1DProfile(GeometryTopBottom.Top, false);
                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>
        /// 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 isBottomAquifer = true, the highest aquifer in the deep cluster of aquifers
        /// - If 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>
        /// <returns>The geometry surface of the selected aquifer (bottom or in-between)</returns>
        private GeometrySurface GeometrySurfaceAquifer(bool isBottomAquifer, bool isTopInBetweenAquifer)
        {
            SoilLayer1D relevantAquiferLayerIn1DProfile = null;
            GeometrySurface relevantSurface = 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;
                }
            }

            var middlePoint = new GeometryPoint();
            middlePoint.X = xDikeTopAtPolder;
            middlePoint.Z = (relevantAquiferLayerIn1DProfile.TopLevel + relevantAquiferLayerIn1DProfile.BottomLevel) / 2 ;
            var relevantAquiferLayerIn2DProfile = profile2D.GetSoilLayer(middlePoint);
            
            return relevantAquiferLayerIn2DProfile.GeometrySurface;
        }

        /// <summary>
        /// order the curves in a geometryloop to make a closed circle
        /// </summary>
        /// <param name="geometryCurves"></param>
        /// <returns></returns>
        private List<GeometryCurve> OrderCurves(GeometryLoop geometryCurves)
        {
            var minX = geometryCurves.GetMinX();
            var minXPoints = geometryCurves.GetPointsAtX(minX);
            var startTopCurves = minXPoints.OrderByDescending(p => p.Z).First();
            var firstCurve = geometryCurves.CurveList.First(c => c.HeadPoint == startTopCurves);
            var orderedCurves = OrderCurves(firstCurve, geometryCurves.CurveList);
            return orderedCurves;
        }

        /// <summary>
        /// create a copy of the curvelist, this will leave the outerloop in its original state
        /// </summary>
        /// <param name="geometrySurface"></param>
        /// <returns></returns>
        private GeometryLoop GetGeometrySurfaceCurvesCopy(GeometrySurface geometrySurface)
        {
            var curvesCopy = new GeometryLoop();
            foreach (var geometryCurve in geometrySurface.OuterLoop.CurveList)
            {
                curvesCopy.CurveList.Add(geometryCurve);
            }
            return curvesCopy;
        }

        /// <summary>
        /// get next curve so that curve.EndPoint == nextCurve.EndPoint
        /// </summary>
        /// <param name="curve"></param>
        /// <param name="boundaryCurves"></param>
        /// <returns></returns>
        private GeometryCurve GetNextCurveReversed(GeometryCurve curve, List<GeometryCurve> boundaryCurves)
        {
            // TODO: Code clone of GeometrySurface.GetNextCurveReversed; Refactoring recommended
            foreach (var geometryCurve in boundaryCurves)
            {
                if (curve.EndPoint == geometryCurve.EndPoint && curve.HeadPoint != geometryCurve.HeadPoint)
                {
                    boundaryCurves.Remove(geometryCurve);
                    return new GeometryCurve(geometryCurve.EndPoint, geometryCurve.HeadPoint);
                }
            }
            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;
        }
    }
}