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

using System;
using System.Collections.Generic;
using System.Linq;
using Deltares.DamEngine.Data.Geometry;
using Deltares.DamEngine.Data.Geotechnics;

namespace Deltares.DamEngine.Calculators.KernelWrappers.Common;

public static class SoilProfile2DHelper
{
    public enum LayerType
    {
        BottomAquiferCluster,
        InBetweenAquiferCluster,
        HighestAquiferCluster,
        LowestLayer
    }

    private enum BoundaryType
    {
        Top,
        Bottom
    }

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

            previousCount = currentCount;
        }

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

            previousCount = currentCount;
            previousCrossSection = crossSection;
        }

        return false;
    }

    private static bool IsAtLeastOneIsolatedInBetweenAquiferPresent(int currentCount, int previousCount, SoilProfile1D currentCrossSection, SoilProfile1D previousCrossSection)
    {
        // If more than one new or one less in-between aquifers are found, at least one is an isolated aquifer.
        if (Math.Abs(currentCount - previousCount) > 1)
        {
            return true;
        }

        if (currentCount <= 0 || currentCrossSection.GetInBetweenAquiferClusters == null)
        {
            return false;
        }

        // If one new in-between aquifer is found, the top layer in the previous cross-section must be an aquifer.
        // If not, it is an isolated aquifer.
        if (currentCount - previousCount == 1 && !previousCrossSection.Layers[0].IsAquifer)
        {
            return true;
        }

        // If one less in-between aquifer is found, the top layer in the current section must be an aquifer.
        // If not, it is an isolated aquifer.
        return previousCount - currentCount == 1 && !currentCrossSection.Layers[0].IsAquifer;
    }

    /// <summary>
    /// Determines all the points of the layer boundary for the specified layer type:
    /// - For LayerType = BottomAquiferCluster, layerBoundaryTop is the top boundary of the bottom aquifer (= waternet line for PL3)
    /// - For LayerType = InBetweenAquiferCluster, layerBoundaryTop and layerBoundaryBottom are both boundaries of the in-between aquifer (= waternet line for PL4) 
    /// - For LayerType = HighestAquiferCluster, layerBoundaryTop is the top of the highest aquifer (= waternet line for PL1 for "Hydrostatic" waternet)
    /// - For LayerType = LowestLayer, layerBoundaryBottom is the bottom boundary of the lowest layer (= waternet line for PL1 for "Full hydrostatic" waternet)
    /// When the aquifer is not continuous, null top and bottom boundaries are returned.
    /// </summary>
    /// <param name="layerType">The type of layer boundary.</param>
    /// <param name="soilProfile">The soil profile 2D.</param>
    /// <param name="layerBoundaryTop">All the points of the layer top boundary.</param>
    /// <param name="layerBoundaryBottom">All the points of the layer bottom boundary.</param>
    /// <param name="indexInBetweenAquifer">In case of several in-between aquifers, the index of the in-between aquifer must be specified.</param>
    /// <returns>All the points of the layer boundary.</returns>
    internal static void DetermineAquiferLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, out Point2D[] layerBoundaryTop, out Point2D[] layerBoundaryBottom, int indexInBetweenAquifer = -1)
    {
        double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile);
        var xCoordinatesAll = new List<double>();
        foreach (double xCoordinate in xCoordinates)
        {
            if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile, indexInBetweenAquifer))
            {
                xCoordinatesAll.Add(xCoordinate - deviationX);
                xCoordinatesAll.Add(xCoordinate + deviationX);
            }
            else
            {
                xCoordinatesAll.Add(xCoordinate);
            }
        }

        var layerBoundaryTopPoints = new List<Point2D>();
        var layerBoundaryBottomPoints = new List<Point2D>();
        double previousAquiferTop = double.NaN;
        double previousAquiferBottom = double.NaN;
        foreach (double xCoordinate in xCoordinatesAll)
        {
            SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate);

            // Determine if the cluster of layers is in range of the previous cluster of layers.
            // If not, return empty coordinates, because the cluster of layers is interrupted.
            double currentAquiferTop = GetLevel(layerType, BoundaryType.Top, crossSection, indexInBetweenAquifer);
            double currentAquiferBottom = GetLevel(layerType, BoundaryType.Bottom, crossSection, indexInBetweenAquifer);

            if (!double.IsNaN(previousAquiferTop) && !double.IsNaN(previousAquiferBottom)
                                                  && !double.IsNaN(currentAquiferTop)
                                                  && !double.IsNaN(currentAquiferBottom)
                                                  && !AreHorizontallyConnected(previousAquiferTop, previousAquiferBottom, currentAquiferTop, currentAquiferBottom))
            {
                layerBoundaryTop = null;
                layerBoundaryBottom = null;
                return;
            }

            if (!double.IsNaN(currentAquiferTop) && !double.IsNaN(currentAquiferBottom))
            {
                previousAquiferTop = currentAquiferTop;
                previousAquiferBottom = currentAquiferBottom;
                layerBoundaryTopPoints.Add(new Point2D(xCoordinate, currentAquiferTop));
                layerBoundaryBottomPoints.Add(new Point2D(xCoordinate, currentAquiferBottom));
            }
        }

        // Perform a short validation that the coordinates are fully defined from the beginning to the end 
        // of the profile. If not, the layer is not continuous.
        if (!IsLayerBoundaryContinuous(layerBoundaryTopPoints, xCoordinates.First(), xCoordinates.Last())
            || !IsLayerBoundaryContinuous(layerBoundaryBottomPoints, xCoordinates.First(), xCoordinates.Last()))
        {
            layerBoundaryTop = null;
            layerBoundaryBottom = null;
            return;
        }

        layerBoundaryTop = layerBoundaryTopPoints.ToArray();
        layerBoundaryBottom = layerBoundaryBottomPoints.ToArray();
    }

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

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

    private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer)
    {
        SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D, indexInBetweenAquifer);
        if (layer != null)
        {
            return boundaryType == BoundaryType.Top ? layer.TopLevel : layer.BottomLevel;
        }

        return double.NaN;
    }

    private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer)
    {
        switch (layerType)
        {
            case LayerType.BottomAquiferCluster:
                return boundaryType == BoundaryType.Top ? soilProfile1D.BottomAquiferLayer : soilProfile1D.DeepestAquiferLayer;
            case LayerType.InBetweenAquiferCluster:
                return boundaryType == BoundaryType.Top ? soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item1 : soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item2;
            case LayerType.HighestAquiferCluster:
                return boundaryType == BoundaryType.Top ? soilProfile1D.GetHighestAquifer() : soilProfile1D.GetLowestLayerOfHighestAquiferCluster();
            case LayerType.LowestLayer:
                return soilProfile1D.Layers.Last();
        }

        return null;
    }

    private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile, int indexInBetweenAquifer)
    {
        SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX);
        SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX);
        if (crossSectionLeft == null || crossSectionRight == null)
        {
            return false;
        }

        double aquiferLeftTop = GetLevel(layerType, BoundaryType.Top, crossSectionLeft, indexInBetweenAquifer);
        double aquiferLefBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionLeft, indexInBetweenAquifer);
        double aquiferRightTop = GetLevel(layerType, BoundaryType.Top, crossSectionRight, indexInBetweenAquifer);
        double aquiferRightBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionRight, indexInBetweenAquifer);
        if (!double.IsNaN(aquiferLeftTop) && !double.IsNaN(aquiferLefBottom) && !double.IsNaN(aquiferRightTop) && !double.IsNaN(aquiferRightBottom))
        {
            return Math.Abs(aquiferLeftTop - aquiferRightTop) > GeometryConstants.Accuracy || Math.Abs(aquiferLefBottom - aquiferRightBottom) > GeometryConstants.Accuracy;
        }

        return false;
    }

    private static bool AreHorizontallyConnected(double leftLayerTop, double leftLayerBottom, double rightLayerTop, double rightLayerBottom)
    {
        // Left soil layer envelopes whole right soil layer
        if (leftLayerBottom <= rightLayerBottom && leftLayerTop >= rightLayerTop)
        {
            return true;
        }

        // Right soil layer envelopes whole left soil layer
        if (rightLayerBottom <= leftLayerBottom && rightLayerTop >= leftLayerTop)
        {
            return true;
        }

        return (rightLayerTop <= leftLayerTop && rightLayerTop >= leftLayerBottom) // Top level lies in-between the left soil layer
               || (rightLayerBottom >= leftLayerBottom && rightLayerBottom <= leftLayerTop); // Bottom level lies in-between the left soil layer
    }
}