// Copyright (C) Stichting Deltares 2025. 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.General;
using Deltares.DamEngine.Data.Geometry;
using Deltares.DamEngine.Data.Geotechnics;
using Deltares.DamEngine.Data.Standard;

namespace Deltares.DamEngine.TestHelpers.Factories;

/// <summary>
/// To indicate the position of the surface line to the soil profile.
/// </summary>
public enum PositionToSoilProfile2D
{
    LeftOfSoilProfile,
    RightOfSoilProfile,
    OnSoilProfile,
    InsideOfSoilProfile
}

public static class FactoryForSoilProfiles
{
    /// <summary>
    /// Create simple soil profile
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateSimpleProfile()
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D();
        layer.TopLevel = 10.0;
        layer.Soil = new Soil("Clay", 16.0, 12.0);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = -20.0;

        return soilProfile;
    }

    /// <summary>
    /// Create two layer soilprofile
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateTwoLayerProfile()
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D();
        layer.TopLevel = 10.0;
        layer.Soil = new Soil("HW-OBO", 12.0, 10.0);
        layer.IsAquifer = false;
        layer.WaterpressureInterpolationModel = WaterpressureInterpolationModel.Automatic;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D();
        layer.TopLevel = 2.0;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        layer.WaterpressureInterpolationModel = WaterpressureInterpolationModel.Hydrostatic;
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = -10.0;

        return soilProfile;
    }

    /// <summary>
    /// Create two layer soil profile with pleistocene sand layer only
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateClaySandProfile(out SoilList soilList, double topLevelClayLayer = 10, double topLevelSandLayer = -5)
    {
        var soilProfile = new SoilProfile1D();

        var soil1 = new Soil("HW-OBO", 12.0, 10.0);
        var soil2 = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        soilList = new SoilList();
        soilList.Add(soil1);
        soilList.Add(soil2);

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = topLevelClayLayer;
        layer.Soil = soil1;
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = topLevelSandLayer;
        layer.Soil = soil2;
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);
        soilProfile.BottomLevel = -10.0;

        return soilProfile;
    }

    /// <summary>
    /// Gammadry and Gammywet are both set to 1 for testing purposes
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateClaySandClaySandProfileForPipingBligh()
    {
        SoilProfile1D soilProfile = CreateClaySandClaySandProfile();
        foreach (SoilLayer1D layer in soilProfile.Layers)
        {
            layer.Soil.BelowPhreaticLevel = 1;
            layer.Soil.AbovePhreaticLevel = 1;
            layer.Soil.DiameterD70 = 0.00018;
        }

        return soilProfile;
    }

    /// <summary>
    /// Creates the clay sand profile for piping bligh.
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateClaySandProfileForPipingBligh(out SoilList soilList)
    {
        SoilProfile1D soilProfile = CreateClaySandProfile(out soilList);
        foreach (SoilLayer1D layer in soilProfile.Layers)
        {
            layer.Soil.BelowPhreaticLevel = 1;
            layer.Soil.AbovePhreaticLevel = 1;
            layer.Soil.DiameterD70 = 0.00018;
        }

        return soilProfile;
    }

    /// <summary>
    /// Create four layer soil profile with pleistocene and intermediate sand layer
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateClaySandClaySandProfile(double topClay1 = 10.0, double topSand1 = 1.0, double topClay2 = -1.5, double topSand2 = -5.0)
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = topClay1;
        layer.Soil = new Soil("HW-OBO", 12.0, 10.0);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = topSand1;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = topClay2;
        layer.Soil = new Soil("HW-DUN", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = topSand2;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = -10.0;
        return soilProfile;
    }

    /// <summary>
    /// Create six layer soil profile with 3 aquifers
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateClaySandClaySandClaySandProfile(double level1 = 10, double level2 = 1, double level3 = -1.5, double level4 = -5, double level5 = -7, double level6 = -5)
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = level1;
        layer.Soil = new Soil("HW-OBO", 12.0, 10.0);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = level2;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = level3;
        layer.Soil = new Soil("HW-DUN", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = level4;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = level5;
        layer.Soil = new Soil("HW-DUN", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = level6;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = level6 - 10;

        return soilProfile;
    }

    /// <summary>
    /// Create 4 layers soil profile with 1 aquifer at the bottom
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D Create4LayersProfileWith1BottomAquifer(double level1 = 0.0, double level2 = -5.0, double level3 = -10.0, double level4 = -15.0, double bottomLevel5 = -20.0, bool isInBetweenAquiferPresent = false)
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile),
            TopLevel = level1,
            IsAquifer = false
        };
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile),
            TopLevel = level2,
            IsAquifer = isInBetweenAquiferPresent
        };
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile),
            TopLevel = level3,
            IsAquifer = false
        };
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile),
            TopLevel = level4,
            IsAquifer = true
        };
        soilProfile.Layers.Add(layer);
        soilProfile.BottomLevel = bottomLevel5;

        return soilProfile;
    }

    /// <summary>
    /// Create four layer soil profile with pleistocene and intermediate sand layer with multiple infiltration layers
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateMultiInfiltrationLayerProfile()
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = 10.0;
        layer.Soil = new Soil("HW-OBO", 12.0, 10.0);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = 1.0;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -1.5;
        layer.Soil = new Soil("HW-DUN", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -2.3;
        layer.Soil = new Soil("HW-OBO", 12.0, 10.0);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -3.2;
        layer.Soil = new Soil("HW-DUN", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -4.1;
        layer.Soil = new Soil("HW-OBO", 12.0, 10.0);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -5.0;
        layer.Soil = new Soil("Alg-zand (0-30)", 22.0, 20.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = -10.0;

        return soilProfile;
    }

    /// <summary>
    /// Create complex soilprofile
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreateComplexProfile()
    {
        var soilProfile = new SoilProfile1D();

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = 10.0;
        layer.Soil = new Soil("HW-OBO", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -0.6;
        layer.Soil = new Soil("HW-DUN", 16.8, 15.8);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -1.8;
        layer.Soil = new Soil("HW-DUOzand", 18.7, 17.7);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -2.8;
        layer.Soil = new Soil("HW-HVN", 11.4, 10.4);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -4.0;
        layer.Soil = new Soil("HW-HVN", 11.40, 10.40);
        layer.IsAquifer = false;
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile)
        };
        layer.TopLevel = -5.0;
        layer.Soil = new Soil("Alg-zand (0-30)", 19.0, 18.0);
        layer.IsAquifer = true;
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = -10.0;
        return soilProfile;
    }
   
    /// <summary>
    /// Create two layer soil profile with pleistocene sand layer (no intermediate layer)
    /// </summary>
    /// <returns>soil profile</returns>
    public static SoilProfile1D CreatePipingSellmeijerProfileWithOneSandLayer(out SoilList soilList)
    {
        var soilProfile = new SoilProfile1D();

        var topMaterial = new Soil("Topmaterial", 14, 11)
        {
            PermeabKx = 0.0003,
            DiameterD70 = Physics.FactorMicroMeterToMeter * 300.0
        };
        var sand = new Soil("Sand", 22.0, 20.0)
        {
            PermeabKx = 0.0001,
            DiameterD70 = Physics.FactorMicroMeterToMeter * 200.0
        };
        soilList = new SoilList();
        soilList.Add(topMaterial);
        soilList.Add(sand);

        var layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile),
            TopLevel = 10.0,
            Soil = topMaterial,
            IsAquifer = false
        };
        soilProfile.Layers.Add(layer);

        layer = new SoilLayer1D
        {
            Name = GetNewUniqueLayerId(soilProfile),
            TopLevel = -2.0,
            Soil = sand,
            IsAquifer = true
        };
        soilProfile.Layers.Add(layer);

        soilProfile.BottomLevel = -10.0;

        return soilProfile;
    }
    
    /// <summary>
    /// Create a soil layer with the given top level and soil name
    /// </summary>
    /// <param name="topLevel">The toplevel.</param>
    /// <param name="soilName">The soil name</param>
    /// <returns></returns>
    public static SoilLayer1D CreateSoilLayer(double topLevel, string soilName)
    {
        var random = new Random(soilName.GetHashCode());

        Array values = Enum.GetValues(typeof(WaterpressureInterpolationModel));
        var model = (WaterpressureInterpolationModel) values.GetValue(random.Next(values.Length));

        return new SoilLayer1D
        {
            TopLevel = topLevel,
            Soil = new Soil(soilName),
            IsAquifer = Convert.ToBoolean(random.Next(0, 2)),
            WaterpressureInterpolationModel = model
        };
    }

    /// <summary>
    /// Create a surface line with the given coordinates
    /// </summary>
    /// <param name="coordinates"></param>
    /// <returns></returns>
    public static SurfaceLine2 CreateSurfaceLine(IEnumerable<Point2D> coordinates)
    {
        var surfaceLine = new SurfaceLine2();
        surfaceLine.Geometry.Points.AddRange(coordinates);

        return surfaceLine;
    }

    /// <summary>
    /// Create a soil profile 2D with two layers
    /// </summary>
    /// <returns>The 2D Soil Profile</returns>
    public static SoilProfile2D CreateSoilProfile2DWithTwoLayers()
    {
        const string layerName = "Layer";
        const string topLayerName = "TopLayer";

        SoilLayer1D soilLayer = CreateSoilLayer(-5, layerName);

        var profile = new SoilProfile1D
        {
            BottomLevel = -10
        };
        profile.Layers.Add(soilLayer);

        SurfaceLine2 surfaceLine = CreateSurfaceLine(new[]
        {
            new Point2D(0, 0),
            new Point2D(5, 10),
            new Point2D(10, 10)
        });

        SoilProfile2D soilProfile2D = CreateSoilProfile2DBasedOnSoilProfile1D(profile, surfaceLine);
        return soilProfile2D;
    }

    /// <summary>
    /// Creates 2D SoilProfile with three layers
    /// </summary>
    /// <returns>The 2D Soil Profile</returns>
    public static SoilProfile2D CreateSoilProfile2DWithThreeLayers()
    {
        const string layer1Name = "Layer1";
        const string layer2Name = "Layer2";
        const string topLayerName = "TopLayer";

        SoilLayer1D soilLayer1 = CreateSoilLayer(-8, layer1Name);
        SoilLayer1D soilLayer2 = CreateSoilLayer(-5, layer2Name);

        var profile = new SoilProfile1D
        {
            BottomLevel = -10
        };
        profile.Layers.Add(soilLayer1);
        profile.Layers.Add(soilLayer2);

        SurfaceLine2 surfaceLine = CreateSurfaceLine(new[]
        {
            new Point2D(0, 0),
            new Point2D(5, 10),
            new Point2D(10, 10)
        });

        SoilProfile2D soilProfile2D = CreateSoilProfile2DBasedOnSoilProfile1D(profile, surfaceLine);
        // Make one boundary not horizontal
        soilProfile2D.Geometry.Points[4].Z = 0.0;
        return soilProfile2D;
    }

    /// <code>
    ///  -50     -30             -10      10             40   50 
    ///   |---------------------------------------------------|  Level 10 m
    ///   |   surface 1                                       | 
    ///   |      |---------------|        |---------------|   |   Level 0 m
    ///   |      |   surface 2   |        |   surface 3   |   |  
    ///   |      |-------------- |        |---------------|   |   Level -10 m
    ///   |   surface 1                                       |
    ///   |---------------------------------------------------|  Level -20 m
    /// </code>
    public static SoilProfile2D CreateSoilProfile2DWithThreeLayersOfWhichTwoAreInnerLoops()
    {
        var soilProfile2D = new SoilProfile2D
        {
            Geometry = new GeometryData
            {
                Left = -50,
                Right = 50,
                Bottom = -20
            }
        };
        var soil1 = new Soil("Soil1", 11, 12);
        var soil2 = new Soil("Soil2", 13, 14);
        var soil3 = new Soil("Soil3", 13.1, 14.1);

        SoilLayer2D soilLayer1 = CreateRectangularSoilLayer2D(10, -20, -50, 50, soilProfile2D, soil1);
        SoilLayer2D soilLayer2 = CreateRectangularSoilLayer2D(0, -10, -30, -10, soilProfile2D, soil2);
        SoilLayer2D soilLayer3 = CreateRectangularSoilLayer2D(0, -10, 10, 40, soilProfile2D, soil3);

        soilProfile2D.Surfaces.Add(soilLayer1);
        soilProfile2D.Surfaces.Add(soilLayer2);
        soilProfile2D.Surfaces.Add(soilLayer3);

        soilProfile2D.Geometry.Loops.Add(soilLayer1.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer1.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer2.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer2.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer3.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer3.GeometrySurface);
        soilProfile2D.Geometry.Surfaces[0].InnerLoops.Add(soilLayer2.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces[0].InnerLoops.Add(soilLayer3.GeometrySurface.OuterLoop);

        return soilProfile2D;
    }

    /// <code>
    ///  -50 -40 -30             0      10       30  40  50 
    ///   |-----------------------------------------------|  Level 10 m
    ///   |   surface 1                                   | 
    ///   |      |---------------|---------------|        |   Level 0 m
    ///   |     |   surface 2    |   surface 3    |       |  
    ///   |   |------------------|------------------|     |   Level -10 m
    ///   |   surface 1                                   |
    ///   |-----------------------------------------------|  Level -20 m
    /// </code>
    public static SoilProfile2D CreateSoilProfile2DWithThreeLayersOfWhichTwoAreConnectedInnerLoops()
    {
        var soilProfile2D = new SoilProfile2D
        {
            Geometry = new GeometryData
            {
                Left = -50,
                Right = 50,
                Bottom = -20
            }
        };
        var soil1 = new Soil("Soil1", 11, 12);
        var soil2 = new Soil("Soil2", 13, 14);
        var soil3 = new Soil("Soil3", 13.1, 14.1);

        SoilLayer2D soilLayer1 = CreateRectangularSoilLayer2D(10, -20, -50, 50, soilProfile2D, soil1);
        SoilLayer2D soilLayer2 = CreateRectangularSoilLayer2D(0, -10, -40, 0, soilProfile2D, soil2);
        SoilLayer2D soilLayer3 = CreateRectangularSoilLayer2D(0, -10, 0, 40, soilProfile2D, soil3);

        soilProfile2D.Geometry.GetPointAtLocation(new Point2D(-40, 0)).X = -30;
        soilProfile2D.Geometry.GetPointAtLocation(new Point2D(40, 0)).X = 30;

        soilProfile2D.Surfaces.Add(soilLayer1);
        soilProfile2D.Surfaces.Add(soilLayer2);
        soilProfile2D.Surfaces.Add(soilLayer3);

        soilProfile2D.Geometry.Loops.Add(soilLayer1.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer1.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer2.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer2.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer3.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer3.GeometrySurface);
        soilProfile2D.Geometry.Surfaces[0].InnerLoops.Add(soilLayer2.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces[0].InnerLoops.Add(soilLayer3.GeometrySurface.OuterLoop);

        return soilProfile2D;
    }

    /// <code>
    /// -50              -20     -10     0              35             60
    /// |-----------------|--------------|-----------------------------|  Level 10 m
    /// |   surface 1     |  surface 2   |        surface 3            |
    /// |-----------------|------|-------|--------------|--------------|  Level 0 m
    /// |       surface 4        |       surface 5      |   surface    |
    /// |                        |                      |      6       |
    /// |------------------------|----------------------|--------------|  Level -15 m
    /// </code>
    public static SoilProfile2D CreateSoilProfile2DWithSixSurfacesFormingTwoLayers(double xLeft = -50)
    {
        var soilProfile2D = new SoilProfile2D
        {
            Geometry = new GeometryData
            {
                Left = xLeft,
                Right = xLeft + 110,
                Bottom = -15
            }
        };
        var soil1 = new Soil("Soil1");
        SoilLayer2D soilLayer1 = CreateRectangularSoilLayer2D(10, 0, xLeft, xLeft + 30, soilProfile2D, soil1);
        var soil2 = new Soil("Soil2");
        SoilLayer2D soilLayer2 = CreatePentagonSoilLayer2D(new Point2D(xLeft + 30, 10), new Point2D(xLeft + 50, 10), new Point2D(xLeft + 50, 0), new Point2D(xLeft + 40, 0), new Point2D(xLeft + 30, 0), soilProfile2D, soil2);
        var soil3 = new Soil("Soil3");
        SoilLayer2D soilLayer3 = CreatePentagonSoilLayer2D(new Point2D(xLeft + 50, 10), new Point2D(xLeft + 110, 10), new Point2D(xLeft + 110, 0), new Point2D(xLeft + 85, 0), new Point2D(xLeft + 50, 0), soilProfile2D, soil3);
        var soil4 = new Soil("Soil4");
        SoilLayer2D soilLayer4 = CreatePentagonSoilLayer2D(new Point2D(xLeft, 0), new Point2D(xLeft + 30, 0), new Point2D(xLeft + 40, 0), new Point2D(xLeft + 40, -15), new Point2D(xLeft, -15), soilProfile2D, soil4);
        var soil5 = new Soil("Soil5");
        SoilLayer2D soilLayer5 = CreatePentagonSoilLayer2D(new Point2D(xLeft + 40, 0), new Point2D(xLeft + 50, 0), new Point2D(xLeft + 85, 0), new Point2D(xLeft + 85, -15), new Point2D(xLeft + 40, -15), soilProfile2D, soil5);
        var soil6 = new Soil("Soil6");
        SoilLayer2D soilLayer6 = CreateRectangularSoilLayer2D(0, -15, xLeft + 85, xLeft + 110, soilProfile2D, soil6, true);

        soilProfile2D.Surfaces.Add(soilLayer1);
        soilProfile2D.Surfaces.Add(soilLayer2);
        soilProfile2D.Surfaces.Add(soilLayer3);
        soilProfile2D.Surfaces.Add(soilLayer4);
        soilProfile2D.Surfaces.Add(soilLayer5);
        soilProfile2D.Surfaces.Add(soilLayer6);

        soilProfile2D.Geometry.Loops.Add(soilLayer1.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer1.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer2.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer2.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer3.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer3.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer4.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer4.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer5.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer5.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer6.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer6.GeometrySurface);

        return soilProfile2D;
    }

    /// <code>
    ///  -50         -20/-19.999/-19.998       -10 
    ///   |-----------------|-|-|---------------|  Level 10 m
    ///   |   surface 1     |2|3|   surface 4   | 
    ///   |-----------------|-|-|---------------|  Level 0 m
    ///   |   surface 5       |6|   surface 7   |
    ///   |                   | |               |
    ///   |-------------------|-|---------------|  Level -15 m
    /// </code>
    public static SoilProfile2D CreateSoilProfile2DWithFourSurfacesAndThreeVerticalSeparationDistantFrom1Mm(bool areConsecutiveSoilsIdentical)
    {
        var soilProfile2D = new SoilProfile2D
        {
            Geometry = new GeometryData
            {
                Left = -50,
                Right = -10,
                Bottom = -15
            }
        };
        var soil1 = new Soil("Soil1", 11, 12);
        Soil soil2 = areConsecutiveSoilsIdentical ? soil1 : new Soil("Soil2", 13, 14);
        Soil soil3 = areConsecutiveSoilsIdentical ? soil1 : new Soil("Soil3", 13.1, 14.1);
        Soil soil4 = areConsecutiveSoilsIdentical ? soil1 : new Soil("Soil4", 13.2, 14.2);
        var soil5 = new Soil("Soil5", 15, 16);
        Soil soil6 = areConsecutiveSoilsIdentical ? soil5 : new Soil("Soil6", 17, 18);
        Soil soil7 = areConsecutiveSoilsIdentical ? soil5 : new Soil("Soil7", 17.1, 18.1);

        SoilLayer2D soilLayer1 = CreateRectangularSoilLayer2D(10, 0, -50, -20, soilProfile2D, soil1);
        SoilLayer2D soilLayer2 = CreateRectangularSoilLayer2D(10, 0, -20, -19.999, soilProfile2D, soil2);
        SoilLayer2D soilLayer3 = CreateRectangularSoilLayer2D(10, 0, -19.999, -19.998, soilProfile2D, soil3);
        SoilLayer2D soilLayer4 = CreateRectangularSoilLayer2D(10, 0, -19.998, -10, soilProfile2D, soil4);
        SoilLayer2D soilLayer5 = CreatePentagonSoilLayer2D(new Point2D(-50, 0), new Point2D(-20, 0), new Point2D(-19.999, 0), new Point2D(-19.999, -15), new Point2D(-50, -15), soilProfile2D, soil5);
        SoilLayer2D soilLayer6 = CreateRectangularSoilLayer2D(0, -15, -19.999, -19.998, soilProfile2D, soil6);
        SoilLayer2D soilLayer7 = CreateRectangularSoilLayer2D(0, -15, -19.998, -10, soilProfile2D, soil7);

        soilProfile2D.Surfaces.Add(soilLayer1);
        soilProfile2D.Surfaces.Add(soilLayer2);
        soilProfile2D.Surfaces.Add(soilLayer3);
        soilProfile2D.Surfaces.Add(soilLayer4);
        soilProfile2D.Surfaces.Add(soilLayer5);
        soilProfile2D.Surfaces.Add(soilLayer6);
        soilProfile2D.Surfaces.Add(soilLayer7);

        soilProfile2D.Geometry.Loops.Add(soilLayer1.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer1.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer2.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer2.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer3.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer3.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer4.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer4.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer5.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer5.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer6.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer6.GeometrySurface);
        soilProfile2D.Geometry.Loops.Add(soilLayer7.GeometrySurface.OuterLoop);
        soilProfile2D.Geometry.Surfaces.Add(soilLayer7.GeometrySurface);

        return soilProfile2D;
    }

    public static SoilProfile1D CreateSoilProfile1DWithTwoClustersOfInBetweenAquifers(out SurfaceLine2 surfaceLine)
    {
        const double topLevel = 1.212;
        surfaceLine = new SurfaceLine2();
        surfaceLine.EnsurePointOfType(0, topLevel, CharacteristicPointType.SurfaceLevelOutside);
        surfaceLine.EnsurePointOfType(25, topLevel, CharacteristicPointType.DikeTopAtRiver);
        surfaceLine.EnsurePointOfType(50, topLevel, CharacteristicPointType.DikeTopAtPolder);
        surfaceLine.EnsurePointOfType(75, topLevel, CharacteristicPointType.DikeToeAtPolder);
        surfaceLine.EnsurePointOfType(100, topLevel, CharacteristicPointType.SurfaceLevelInside);

        var soilProfile1D = new SoilProfile1D
        {
            BottomLevel = -30.0
        };

        var layer = new SoilLayer1D
        {
            Name = "L1a",
            TopLevel = topLevel,
            IsAquifer = true // aquifer at top
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L1b",
            TopLevel = 0,
            IsAquifer = false
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L2a",
            TopLevel = -2.111, // top of highest in-between aquifer 
            IsAquifer = true
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L2b",
            TopLevel = -2.151,
            IsAquifer = true
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L3",
            TopLevel = -3.373,
            IsAquifer = false
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L4",
            TopLevel = -4.151, // top of lowest in-between aquifer 
            IsAquifer = true
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L5",
            TopLevel = -5.373,
            IsAquifer = false
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L6a",
            TopLevel = -6.111, // top of bottom aquifer
            IsAquifer = true
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L6b",
            TopLevel = -6.151,
            IsAquifer = true
        };
        soilProfile1D.Layers.Add(layer);
        layer = new SoilLayer1D
        {
            Name = "L7",
            TopLevel = -7.373,
            IsAquifer = false
        };
        soilProfile1D.Layers.Add(layer);

        return soilProfile1D;
    }
    
    /// <code>
    ///                /------\
    ///               /        \
    /// _____________/          \______________________________  Level 10 m
    ///                Aquitard    
    /// _______________________________________________________  Level 6 m                                
    ///             In between aquifer  
    /// _______________________________________________________  Level 2 m          
    ///               Aquitard
    /// _______________________________________________________  Level -2 m
    ///            In between aquifer
    ///              ______________                              Level -3 m
    /// ____________/              \___________________________  Level -4 m 
    ///                Aquitard
    /// _______________________________________________________  Level -6 m
    ///               Bottom aquifer                      
    /// _______________________________________________________  Level -10 m
    /// </code>
    public static SoilProfile2D CreateSoilProfile2DWithNonHorizontalBottomAquitard(out SurfaceLine2 surfaceLine)
    {
        const double almostEqual = 1e-09;
        
        surfaceLine = FactoryForSurfaceLines.CreateSurfaceLineDike(10);
        SoilProfile1D soilProfile1D = CreateClaySandClaySandClaySandProfile(10, 6, 2, -2, -4, -6);
        soilProfile1D.BottomLevel = -10;
        SoilProfile2D soilProfile2D = FactoryForSoilProfiles.CreateSoilProfile2DBasedOnSoilProfile1D(soilProfile1D, surfaceLine);
        Point2D point1 = soilProfile2D.Geometry.Points.First(p => Math.Abs(p.X - 34.5) < almostEqual && Math.Abs(p.Z + 4) < almostEqual);
        point1.Z = -3;
        Point2D point2 = soilProfile2D.Geometry.Points.First(p => Math.Abs(p.X - 40.5) < almostEqual && Math.Abs(p.Z + 4) < almostEqual);
        point2.Z = -3;
        return soilProfile2D;
    }

    public static SoilLayer2D CreateRectangularSoilLayer2D(double topCoord, double bottomCoord, double leftCoord,
                                                           double rightCoord, SoilProfile2D soilProfile2D, Soil soil = null, bool isAquifer = false)
    {
        var topLeftPoint = new Point2D(leftCoord, topCoord);
        var topRightPoint = new Point2D(rightCoord, topCoord);
        var bottomRightPoint = new Point2D(rightCoord, bottomCoord);
        var bottomLeftPoint = new Point2D(leftCoord, bottomCoord);
        var curve1 = new GeometryCurve(topLeftPoint, topRightPoint);
        var curve2 = new GeometryCurve(topRightPoint, bottomRightPoint);
        var curve3 = new GeometryCurve(bottomRightPoint, bottomLeftPoint);
        var curve4 = new GeometryCurve(bottomLeftPoint, topLeftPoint);

        if (soilProfile2D != null)
        {
            AddPointIfNotYetPresentInGeometry(soilProfile2D, topLeftPoint);
            AddPointIfNotYetPresentInGeometry(soilProfile2D, topRightPoint);
            AddPointIfNotYetPresentInGeometry(soilProfile2D, bottomRightPoint);
            AddPointIfNotYetPresentInGeometry(soilProfile2D, bottomLeftPoint);
            AddCurveIfNotYetPresentInGeometry(soilProfile2D, curve1);
            AddCurveIfNotYetPresentInGeometry(soilProfile2D, curve2);
            AddCurveIfNotYetPresentInGeometry(soilProfile2D, curve3);
            AddCurveIfNotYetPresentInGeometry(soilProfile2D, curve4);
        }

        var soilLayer2D = new SoilLayer2D
        {
            GeometrySurface = new GeometrySurface
            {
                OuterLoop = new GeometryLoop
                {
                    CurveList =
                    {
                        curve1,
                        curve2,
                        curve3,
                        curve4
                    }
                }
            },
            Soil = soil,
            IsAquifer = isAquifer
        };
        return soilLayer2D;
    }

    public static SoilLayer2D CreateTriangularSoilLayer2D(Point2D point1, Point2D point2, Point2D point3, SoilProfile2D soilProfile2D, Soil soil = null, bool isAquifer = false)
    {
        return CreatePolygoneSoilLayer2D([
            ..new[]
            {
                point1,
                point2,
                point3
            }
        ], soil, soilProfile2D, isAquifer);
    }

    public static SoilLayer2D CreateQuadrilateralSoilLayer2D(Point2D point1, Point2D point2, Point2D point3, Point2D point4, SoilProfile2D soilProfile2D, Soil soil = null, bool isAquifer = false)
    {
        return CreatePolygoneSoilLayer2D([
            ..new[]
            {
                point1,
                point2,
                point3,
                point4
            }
        ], soil, soilProfile2D, isAquifer);
    }

    public static SoilLayer2D CreatePentagonSoilLayer2D(Point2D point1, Point2D point2, Point2D point3, Point2D point4,
                                                        Point2D point5, SoilProfile2D soilProfile2D, Soil soil = null, bool isAquifer = false)
    {
        return CreatePolygoneSoilLayer2D([
            ..new[]
            {
                point1,
                point2,
                point3,
                point4,
                point5
            }
        ], soil, soilProfile2D, isAquifer);
    }

    public static SoilLayer2D CreateHexagonSoilLayer2D(Point2D point1, Point2D point2, Point2D point3, Point2D point4, Point2D point5, Point2D point6, SoilProfile2D soilProfile2D, Soil soil = null, bool isAquifer = false)
    {
        return CreatePolygoneSoilLayer2D([
            ..new[]
            {
                point1,
                point2,
                point3,
                point4,
                point5,
                point6
            }
        ], soil, soilProfile2D, isAquifer);
    }

    public static SoilLayer2D CreateHeptagonSoilLayer2D(Point2D point1, Point2D point2, Point2D point3, Point2D point4, Point2D point5, Point2D point6, Point2D point7, SoilProfile2D soilProfile2D, Soil soil = null, bool isAquifer = false)
    {
        return CreatePolygoneSoilLayer2D([
            ..new[]
            {
                point1,
                point2,
                point3,
                point4,
                point5,
                point6,
                point7
            }
        ], soil, soilProfile2D, isAquifer);
    }

    public static SoilLayer2D CreatePolygoneSoilLayer2D(List<Point2D> points, Soil soil, SoilProfile2D soilProfile2D, bool isAquifer = false)
    {
        var soilLayer2D = new SoilLayer2D
        {
            GeometrySurface = new GeometrySurface
            {
                OuterLoop = new GeometryLoop()
            },
            Soil = soil,
            IsAquifer = isAquifer
        };
        var curves = new List<GeometryCurve>();
        for (var i = 0; i < points.Count - 1; i++)
        {
            var curve = new GeometryCurve(points[i], points[i + 1]);
            curves.Add(curve);
            soilLayer2D.GeometrySurface.OuterLoop.CurveList.Add(curve);
        }

        var endCurve = new GeometryCurve(points[^1], points[0]);
        curves.Add(endCurve);
        soilLayer2D.GeometrySurface.OuterLoop.CurveList.Add(endCurve);
        if (soilProfile2D != null)
        {
            foreach (Point2D point in points)
            {
                AddPointIfNotYetPresentInGeometry(soilProfile2D, point);
            }

            foreach (GeometryCurve curve in curves)
            {
                AddCurveIfNotYetPresentInGeometry(soilProfile2D, curve);
            }
        }

        return soilLayer2D;
    }

    /// <summary>
    /// Creates a surface line for a soil profile 2D at the requested position.
    /// </summary>
    /// <param name="soilProfile2D"></param>
    /// <param name="positionToSoilProfile2D"></param>
    /// <returns></returns>
    public static SurfaceLine2 CreateSurfaceLineForSoilProfile2D(SoilProfile2D soilProfile2D, PositionToSoilProfile2D positionToSoilProfile2D)
    {
        Point2D geometryPoint = soilProfile2D.Geometry.SurfaceLine.Points.First();
        var leftPoint = new Point2D(geometryPoint.X, geometryPoint.Z);
        geometryPoint = soilProfile2D.Geometry.SurfaceLine.Points.Last();
        var rightPoint = new Point2D(geometryPoint.X, geometryPoint.Z);
        var middlePoint = new Point2D((leftPoint.X + rightPoint.X) / 2, (leftPoint.Z + rightPoint.Z) / 2);
        switch (positionToSoilProfile2D)
        {
            case PositionToSoilProfile2D.LeftOfSoilProfile:
                leftPoint.X -= 1;
                break;
            case PositionToSoilProfile2D.RightOfSoilProfile:
                rightPoint.X += 1;
                break;
            case PositionToSoilProfile2D.OnSoilProfile:
                break;
            case PositionToSoilProfile2D.InsideOfSoilProfile:
                leftPoint.X += 1;
                rightPoint.X -= 1;
                break;
        }

        SurfaceLine2 surfaceLine = CreateSurfaceLine(new[]
        {
            leftPoint,
            middlePoint,
            rightPoint
        });
        surfaceLine.CharacteristicPoints.Add(new CharacteristicPoint(surfaceLine.CharacteristicPoints, leftPoint));
        surfaceLine.CharacteristicPoints.Annotate(0, CharacteristicPointType.SurfaceLevelOutside);
        surfaceLine.CharacteristicPoints.Add(new CharacteristicPoint(surfaceLine.CharacteristicPoints, rightPoint));
        surfaceLine.CharacteristicPoints.Annotate(1, CharacteristicPointType.SurfaceLevelInside);
        return surfaceLine;
    }

    public static SoilProfile2D CreateSoilProfile2DBasedOnSoilProfile1D(SoilProfile1D soilProfile1D, SurfaceLine2 surfaceLine)
    {
        var soilSurfaceProfile = new SoilSurfaceProfile
        {
            SoilProfile = soilProfile1D,
            SurfaceLine2 = surfaceLine,
            DikeEmbankmentMaterial = soilProfile1D.Layers.OrderBy(l => l.TopLevel).First().Soil,
            Name = "Test"
        };
        return soilSurfaceProfile.ConvertToSoilProfile2D();
    }

    private static string GetNewUniqueLayerId(SoilProfile1D soilProfile1D)
    {
        var num = 0;
        string id;
        do
        {
            id = "L" + num++;
        } while (soilProfile1D.GetLayerWithName(id) != null);

        return id;
    }

    private static void AddPointIfNotYetPresentInGeometry(SoilProfile2D soilProfile, Point2D point)
    {
        if (!soilProfile.Geometry.Points.Any(p => p.X.IsNearEqual(point.X) && p.Z.IsNearEqual(point.Z)))
        {
            soilProfile.Geometry.Points.Add(point);
        }
    }

    private static void AddCurveIfNotYetPresentInGeometry(SoilProfile2D soilProfile, GeometryCurve curve)
    {
        if (!soilProfile.Geometry.Curves.Any(c => (c.HeadPoint.X.IsNearEqual(curve.HeadPoint.X) && c.HeadPoint.Z.IsNearEqual(curve.HeadPoint.Z)
                                                                                                && c.EndPoint.X.IsNearEqual(curve.EndPoint.X) && c.EndPoint.Z.IsNearEqual(curve.EndPoint.Z)) ||
                                                  (c.EndPoint.X.IsNearEqual(curve.HeadPoint.X) && c.EndPoint.Z.IsNearEqual(curve.HeadPoint.Z) &&
                                                   c.HeadPoint.X.IsNearEqual(curve.EndPoint.X) && c.HeadPoint.Z.IsNearEqual(curve.EndPoint.Z))))
        {
            soilProfile.Geometry.Curves.Add(curve);
        }
    }
}