using System;
using System.Collections.Generic;

// =======================================================================================================================
// Class name:
//
// Description:
//
// Copyright (c) 2009 Deltares
//
// Date        ID   Modification
// 2009-04-07       Created
// =======================================================================================================================

namespace Deltares.Stability.Calculation.Inner
{
    public class TZones
    {
        // zone 1 and 2 data
        private double FDikeTableHeight = 0;
        private double FDikeTableWidth = 3;
        private TInterfaceData FInterfaceData = null;
        private TLayerRecord[] FLayers;
        private List<TRestProfileRec> FRestProfile = new List<TRestProfileRec>();
        private List<TRestProfileRec> FSafeProfile = new List<TRestProfileRec>();
        private double FSafetyZone1a = 0;
        private double FSafetyZone1b = 0;
        private double FSafetyZone2a = 0;
        private double FSafetyZone2b = 0;
        private double FSafetyZone3a = 0;
        private double FSafetyZone3b = 0;
        private TStabilitySideSet FStabilitySide;
        private bool FStandardOvertopping = false;
        private TPoints[] FSurface;
        // zone 3 data
        private double FXBotLeft = 0;
        private double FXBotRight = 0;
        private double FXMHWInfluence = 0;
        private double FXRoadLeft = 0;
        private double FXRoadRight = 0;
        private double FXStart = 0;
        private double FYBotLeft = 0;
        private double FYBotRight = 0;
        private double FYMHWInfluence = 0;
        private double FYRoadLeft = 0;
        private double FYRoadRight = 0;
        private double XSafeprofileStart;
        private double XStartDOV;
        private double YStartDOV;
        private bool isDOVCalculation;
        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-07       Created
        // =======================================================================================================================
        //Constructor  Create()
        public TZones() : base()
        {
            // TODO: Add any constructor code here
        }

        public TInterfaceData InterfaceData
        {
            get
            {
                return FInterfaceData;
            }
            set
            {
                FInterfaceData = value;
            }
        }

        public List<TRestProfileRec> RestProfile
        {
            get
            {
                return FRestProfile;
            }
            set
            {
                FRestProfile = value;
            }
        }

        public List<TRestProfileRec> SafeProfile
        {
            get
            {
                return FSafeProfile;
            }
            set
            {
                FSafeProfile = value;
            }
        }

        public bool CreateRestProfile()
        {
            InitializeLocalZoneData();

            if (isDOVCalculation)
            {
                if (CreateDOVRestprofile())
                {
                    CreateSafeProfile();
                    return true;
                }
            }
            bool result;
            const double CCoorAccuracy = 0.001;
            int Number;
            double xcoor;
            double ZCoor;
            double ZBot;
            double xLimit;
            int Lsign;
            bool LIsInAir = false;
            // if point in the air then a vertical line is followed
            int LoopCount;
            int LoopMax;
            double ClaySlope;
            double SandSlope;
            bool LFound = false;
            ClaySlope = 2;
            SandSlope = 4;
            LoopMax = 100;
            result = false;
            ZBot = DetermineLowestCoordinateFromLayers();
            Number = 0;
            if ((FStabilitySide == TStabilitySideSet.LeftSide))
            {
                Lsign = -1;
            }
            else
            {
                Lsign = 1;
            }
            if ((Lsign > 0))
            {
                xLimit = DetermineMostRightCoordinateFromLayers();
            }
            else
            {
                xLimit = DetermineMostLeftCoordinateFromLayers();
            }
            // Find first point (x at diketableheigth)
            // this point must be in a soillayer
            xcoor = FXStart;
            ZCoor = FDikeTableHeight;
            AddRestprofilePoint(ref Number, xcoor, ZCoor);
            xcoor = xcoor + Lsign*FDikeTableWidth;
            AddRestprofilePoint(ref Number, xcoor, ZCoor);
            if ((IsPointInAir(xcoor, ZCoor)))
            {
                ZCoor = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, xcoor);
                AddRestprofilePoint(ref Number, xcoor, ZCoor);
            }
            LoopCount = 0;
            do
            {
                LoopCount ++;
                FindIntersectionPoint(Lsign, ClaySlope, xLimit, ZBot, ref xcoor, ref ZCoor, ref LFound, ref LIsInAir);
                if (LFound)
                {
                    // first see if the restprofile goes to clay
                    AddRestprofilePoint(ref Number, xcoor, ZCoor);
                }
                if (!LFound && !LIsInAir)
                {
                    // If not so the restprofile passes sand
                    FindIntersectionPoint(Lsign, SandSlope, xLimit, ZBot, ref xcoor, ref ZCoor, ref LFound, ref LIsInAir);
                    if (LFound)
                    {
                        AddRestprofilePoint(ref Number, xcoor, ZCoor);
                    }
                }
                if (LIsInAir)
                {
                    // the dike is followed
                    MStabDatafunctions.GetnextpointOnSurface(FSurface, Lsign, ref xcoor, ref ZCoor);
                    AddRestprofilePoint(ref Number, xcoor, ZCoor);
                }
                if (!LFound && !LIsInAir)
                {
                    // Under a Clay slope sand is found while under a sand slope clay is found
                    // follow the boundary between those layers to the next point
                }
            } while (!((Math.Abs(xcoor - xLimit) < CCoorAccuracy) || (Math.Abs(ZCoor - ZBot) < CCoorAccuracy) || (LoopCount >= LoopMax)));
            if ((LoopCount < LoopMax))
            {
                result = true;
            }
            // Bepaal ook de grensen van zone 3
            if ((Math.Abs(FXRoadLeft - FXRoadRight) > 0.1))
            {
                // er is een zone 3 aanwezig
                // zet de zone grensen eventuel goed
                if (FXRoadLeft > FXRoadRight)
                {
                    MStabDatafunctions.SwapDoubles(ref FXRoadLeft, ref FXRoadRight);
                }
                // Bereken de zone grensen
                FYRoadLeft = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, FXRoadLeft);
                FYRoadRight = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, FXRoadRight);
                FYBotLeft = ZBot;
                FYBotRight = ZBot;
                // Helling 1:1
                FXBotLeft = FXRoadLeft - Math.Abs(FYRoadLeft - FYBotLeft);
                FXBotRight = FXRoadRight - Math.Abs(FYRoadRight - FYBotLeft);
            }
            else
            {
                // geen zone 3, zet alle waarden op 0
                FYRoadLeft = 0;
                FYRoadRight = 0;
                FXBotLeft = 0;
                FXBotRight = 0;
                FYBotLeft = 0;
                FYBotRight = 0;
            }
            return result;
        }

        public void GetZoneNr(bool AIsBishopCalculation, double AXMid, double AZMid, double ARad, double AFSafety,
                              double AXSurface, ref double Angle, ref int AZoneNr, ref int AExtraZone, double AZExit,
                              ref bool APlaneIntesectsRestProfile)
        {
            double LZSurface = 0;
            double LXRestProfile = 0;
            double LZRestProfile = 0;
            double LZSafeLine = 0;
            int LSign = 0;
            double LAngle1 = 0;
            double angle2 = 0;
            double LMinRotation = 0;
            bool LContinue = true;
            if (FRestProfile.Count == 0)
            {
                CreateRestProfile();
            }
            APlaneIntesectsRestProfile = CirkelIntersectsRestProfile(AXMid, AZMid, ARad, ref LXRestProfile, ref LZRestProfile);
            if (isDOVCalculation)
            {
                if (FSafeProfile.Count == 0)
                {
                    CreateSafeProfile();
                }
                if (DetermineZSafeprofileAtX(LXRestProfile, ref LZSafeLine))
                {
                    if (LZSafeLine > LZRestProfile)
                    {
                        AZoneNr = 3; // zone 2
                        LContinue = false;
                    }
                }
            }
            if (LContinue)
            {
                if (APlaneIntesectsRestProfile)
                {
                    LMinRotation = Constants.CMinRotation*10;
                }
                else
                {
                    LMinRotation = 0;
                }
                Angle = 0;
                AZoneNr = 1;
                AExtraZone = 0;
                if ((FStabilitySide == TStabilitySideSet.LeftSide))
                {
                    LSign = -1;
                }
                else
                {
                    LSign = 1;
                }
                // bepaal de Y van het intree punt in het dijk lichaam
                if ((Math.Abs(AXSurface - AXMid) <= ARad + Constants.CGeoAccu))
                {
                    if ((Math.Abs(Math.Abs(AXSurface - AXMid) - ARad) < Constants.CGeoAccu))
                    {
                        LZSurface = AZMid;
                    }
                    else
                    {
                        // de cirkel moet in ieder geval het maaiveld snijden
                        LZSurface = AZMid - Math.Sqrt(ARad*ARad - (AXSurface - AXMid)*(AXSurface - AXMid));
                    }
                    // Bepaal of xsurf en YSurf op rest profiel liggen
                    if (IsPointOnRestProfile(AXSurface, LZSurface))
                    {
                        AZoneNr = 1;
                    }
                    else
                    {
                        // In principe een "zone twee" cirkel, maar kijken of het restprofile wordt gesneden
                        // in geval van mlift kan nog niet naar de veilighid gekeken worden
                        // omdat meerdere glijvlakken hetzelfde middelpunt en entree punt hebben
                        AZoneNr = 3;
                        Angle = LMinRotation;
                        if ((AFSafety <= FSafetyZone1a) || (AFSafety <= FSafetyZone1b))
                        {
                            // (*  or (aFSafety <= FSafetyZone1b) or (not AIsBishopCalculation) then  *)
                            if (APlaneIntesectsRestProfile)
                            {
                                if ((LZRestProfile < AZExit))
                                {
                                    // in such a case rotated mass will always be stable
                                    // meaning zone nr stays 3
                                    Angle = 0;
                                }
                                else
                                {
                                    // Calculate the rotation angle
                                    LAngle1 = Math.Atan2((LZSurface - AZMid), (AXSurface - AXMid));
                                    angle2 = Math.Atan2((LZRestProfile - AZMid), (LXRestProfile - AXMid));
                                    Angle = LAngle1 - angle2;
                                    if (Math.Abs(Angle) < LMinRotation)
                                    {
                                        Angle = Math.Sign(Angle)*LMinRotation;
                                    }
                                }
                            }
                        }
                    }
                    // kijk of de cirkel ook een zone 3 cirkel is (def
                    AExtraZone = AZoneNr;
                    if ((CirkelIntersectsZone3(AXMid, AZMid, ARad)))
                    {
                        AExtraZone = 5;
                    }
                    // Kijk cirkel invloed van MHW heeft
                    if ((LSign*AXSurface > LSign*FXMHWInfluence) && ((AZMid - ARad) > FYMHWInfluence))
                    {
                        AZoneNr = AZoneNr + 1;
                        AExtraZone = AExtraZone + 1;
                    }
                }
            }
        }

        /// <summary>
        /// Determines the zone data for spencer calculation
        /// </summary>
        /// <param name="ACoorArray"></param>
        /// <param name="AYMin"></param>
        /// <returns></returns>
        public int GetSpencerZoneNr(List<TPoints> ACoorArray)
        {
            int i;
            int result;
            double LXS;
            double LYS;
            double YMin = 0;

            if (FRestProfile.Count == 0)
            {
                CreateRestProfile();
            }

            result = 3;

            if ((FStabilitySide == TStabilitySideSet.LeftSide))
            {
                LXS = ACoorArray[ACoorArray.Count - 1].XCoor;
                LYS = ACoorArray[ACoorArray.Count - 1].ZCoor;
            }
            else
            {
                LXS = ACoorArray[0].XCoor;
                LYS = ACoorArray[0].ZCoor;
            }

            if (PointInRestProfile(LXS, LYS))
            {
                result = 1;
            }
            else
            {
                for (i = 0; i < ACoorArray.Count - 1; i++)
                {
                    {
                        if (LineIntersectsRestProfile(ACoorArray[i], ACoorArray[i + 1], ref LXS, ref LYS))
                        {
                            if (DetermineZSafeprofileAtX(LXS, ref YMin))
                            {
                                if (LYS > YMin)
                                {
                                    result = 1;
                                    break;
                                }
                            }
                        }
                    }
                }
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-09       Created
        // =======================================================================================================================
        private void InitializeLocalZoneData()
        {
            // create short cutsfor use in this unit
            FDikeTableHeight = FInterfaceData.ZoneData.DikeTableHeight;
            FDikeTableWidth = FInterfaceData.ZoneData.DikeTableWidth;
            FXStart = FInterfaceData.ZoneData.XStart;
            FXMHWInfluence = FInterfaceData.ZoneData.XMHWInfluence;
            FYMHWInfluence = FInterfaceData.ZoneData.YMHWInfluence;
            FSafetyZone1a = FInterfaceData.ZoneData.SafetyZone1a;
            FSafetyZone1b = FInterfaceData.ZoneData.SafetyZone1b;
            FSafetyZone2a = FInterfaceData.ZoneData.SafetyZone2a;
            FSafetyZone2b = FInterfaceData.ZoneData.SafetyZone2b;
            // zone 3 data
            FXRoadLeft = FInterfaceData.ZoneData.XRoadLeft;
            FXRoadRight = FInterfaceData.ZoneData.XRoadRight;
            FYRoadLeft = FInterfaceData.ZoneData.YRoadLeft;
            FYRoadRight = FInterfaceData.ZoneData.YRoadRight;
            FXBotLeft = FInterfaceData.ZoneData.XBotLeft;
            FXBotRight = FInterfaceData.ZoneData.XBotRight;
            FYBotLeft = FInterfaceData.ZoneData.YBotLeft;
            FYBotRight = FInterfaceData.ZoneData.YBotRight;
            FSafetyZone3a = FInterfaceData.ZoneData.SafetyZone3a;
            FSafetyZone3b = FInterfaceData.ZoneData.SafetyZone3b;

            // DOV data
            isDOVCalculation = FInterfaceData.ZoneData.IsDOVCalculation;
            XStartDOV = FInterfaceData.ZoneData.XStartDOV;
            YStartDOV = FInterfaceData.ZoneData.YStartDOV;
            XSafeprofileStart = FInterfaceData.ZoneData.XSafeprofileStart;

            // standard data
            FStabilitySide = FInterfaceData.ZoneData.StabilitySide;
            FStandardOvertopping = FInterfaceData.ZoneData.StandardOvertopping;
            FLayers = FInterfaceData.Layers;
            FSurface = FInterfaceData.SurfaceLine;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-09       Created
        // =======================================================================================================================
        private double DetermineLowestCoordinateFromLayers()
        {
            double result;
            int i;
            int LNDim;
            result = 0;
            LNDim = FLayers.Length;
            if ((LNDim > 0))
            {
                result = FLayers[0].MinZ;
                for (i = 1; i < LNDim; i ++)
                {
                    result = Math.Min(result, FLayers[i].MinZ);
                }
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-15       Created
        // =======================================================================================================================
        private double DetermineMostRightCoordinateFromLayers()
        {
            double result;
            int i;
            int LNDim;
            result = 100000;
            LNDim = FLayers.Length;
            if ((LNDim > 0))
            {
                result = FLayers[0].MaxX;
                for (i = 1; i < LNDim; i ++)
                {
                    result = Math.Max(result, FLayers[i].MaxX);
                }
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-15       Created
        // =======================================================================================================================
        private double DetermineMostLeftCoordinateFromLayers()
        {
            double result;
            int i;
            int LNDim;
            result = -1000000;
            LNDim = FLayers.Length;
            if ((LNDim > 0))
            {
                result = FLayers[0].MinX;
                for (i = 1; i < LNDim; i ++)
                {
                    result = Math.Min(result, FLayers[i].MinX);
                }
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-15       Created
        // =======================================================================================================================
        private void AddRestprofilePoint(ref int number, double xcoor, double ycoor)
        {
            const double CCoorAccuracy = 0.001;
            double LAccu;
            LAccu = CCoorAccuracy;
            if ((number == 0))
            {
                number ++;
                var RestRec = new TRestProfileRec();
                RestRec.XCoor = xcoor;
                RestRec.YCoor = ycoor;
                FRestProfile.Add(RestRec);
            }
            else
            {
                if ((Math.Abs(xcoor - RestProfile[number - 1].XCoor) > LAccu) ||
                    (Math.Abs(ycoor - RestProfile[number - 1].YCoor) > LAccu))
                {
                    // Niet hetzelfde punt nogmaals toevoegen
                    number ++;
                    var RestRec = new TRestProfileRec();
                    RestRec.XCoor = xcoor;
                    RestRec.YCoor = ycoor;
                    FRestProfile.Add(RestRec);
                }
            }
        }

        private void AddSafeprofilePoint(ref int number, double xcoor, double ycoor)
        {
            const double CCoorAccuracy = 0.001;
            double LAccu;
            LAccu = CCoorAccuracy;
            if ((number == 0))
            {
                number ++;
                var RestRec = new TRestProfileRec();
                RestRec.XCoor = xcoor;
                RestRec.YCoor = ycoor;
                FSafeProfile.Add(RestRec);
            }
            else
            {
                if ((Math.Abs(xcoor - FSafeProfile[number - 1].XCoor) > LAccu) ||
                    (Math.Abs(ycoor - FSafeProfile[number - 1].YCoor) > LAccu))
                {
                    // Niet hetzelfde punt nogmaals toevoegen
                    number ++;
                    var RestRec = new TRestProfileRec();
                    RestRec.XCoor = xcoor;
                    RestRec.YCoor = ycoor;
                    FSafeProfile.Add(RestRec);
                }
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-15       Created
        // =======================================================================================================================
        private bool IsPointInAir(double XCoor, double ZCoor)
        {
            bool result;
            result = (MStabDatafunctions.FindLayerNumberFromPoint(FLayers, XCoor, ZCoor) == -1);
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-07-08       Created
        // =======================================================================================================================
        private void RemoveExcistingPoint(double XCoor, double ZCoor, List<double> LXIntersect, List<double> LZIntersect)
        {
            const double cCooraccuracy = 0.001;

            for (int i = LXIntersect.Count - 1; i >= 0; i--)
            {
                if (((Math.Abs(LXIntersect[i] - XCoor) < cCooraccuracy) &&
                     (Math.Abs(LZIntersect[i] - ZCoor) < cCooraccuracy)))
                {
                    LXIntersect.RemoveAt(i);
                    LZIntersect.RemoveAt(i);
                }
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-16       Created
        // =======================================================================================================================
        private void FindIntersectionPoint(int ASign, double ASlope, double AXLimit, double AZBot, ref double XCoor,
                                           ref double ZCoor, ref bool Found, ref bool IsInAir)
        {
            const double CMinSlopeDiff = 0.01;
            double LXEnd = 0;
            double LSXCoor = 0;
            double LSZCoor = 0;
            int LNIntersect = 0;
            var LXIntersect = new List<double>();
            var LZIntersect = new List<double>();
            int LayerNumber = 0;
            Found = false;
            IsInAir = true;
            // Determine a fictive end point
            LXEnd = XCoor + ASign*Math.Abs(ZCoor - AZBot)*ASlope;
            // Find the firt intersection point with the layers
            MStabDatafunctions.IntersectLineSoileAreas(FLayers, XCoor, ZCoor, LXEnd, AZBot, LXIntersect, LZIntersect);
            RemoveExcistingPoint(XCoor, ZCoor, LXIntersect, LZIntersect);
            if ((LXIntersect.Count > 0))
            {
                // The intersection points are sorted from low x to high x
                if ((ASign > 0))
                {
                    LSXCoor = LXIntersect[0];
                    LSZCoor = LZIntersect[0];
                }
                else
                {
                    LSXCoor = LXIntersect[LNIntersect - 1];
                    LSZCoor = LZIntersect[LNIntersect - 1];
                }
                // Find material from thet layer
                LayerNumber = MStabDatafunctions.FindLayerNumberFromPoint(FLayers, (XCoor + LSXCoor)/2,
                                                                          (ZCoor + LSZCoor)/2);
                IsInAir = (LayerNumber == -1);
                if (!IsInAir)
                {
                    Found =
                        Math.Abs(InterfaceData.StabSoilData[FLayers[LayerNumber].MaterialNumber].RestSlope - ASlope) <
                        CMinSlopeDiff;
                }
                if (Found)
                {
                    XCoor = LSXCoor;
                    ZCoor = LSZCoor;
                }
            }
        }

        private void CreateSafeProfile()
        {
            int i;
            int LNumber = 0;
            double LLevel;
            double pointx;
            double pointz;
            double zcoor;

            LLevel = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, XSafeprofileStart);
            AddSafeprofilePoint(ref LNumber, XSafeprofileStart, LLevel);
            if (FStabilitySide == TStabilitySideSet.RightSide)
            {
                int maxpoint = FSurface.GetUpperBound(0);
                int minpoint = FSurface.GetLowerBound(0);
                for (i = maxpoint; i >= minpoint; i--)
                {
                    pointx = FSurface[i].XCoor;
                    pointz = FSurface[i].ZCoor;
                    if (pointx < XSafeprofileStart)
                    {
                        if (pointz > LLevel)
                        {
                            zcoor = 0.5*(LLevel + pointz);
                            AddSafeprofilePoint(ref LNumber, pointx, zcoor);
                        }
                        else
                        {
                            AddSafeprofilePoint(ref LNumber, pointx, pointz);
                        }
                    }
                }
            }
            else
            {
                for (i = FSurface.GetLowerBound(0) + 1; i <= FSurface.GetUpperBound(0); i++)
                {
                    pointx = FSurface[i].XCoor;
                    if (pointx > XSafeprofileStart)
                    {
                        pointz = FSurface[i].ZCoor;
                        if (pointz > LLevel)
                        {
                            zcoor = 0.5*(LLevel + pointz);
                            AddSafeprofilePoint(ref LNumber, pointx, zcoor);
                        }
                        else
                        {
                            AddSafeprofilePoint(ref LNumber, pointx, pointz);
                        }
                    }
                }
            }
        }

        private void addPointsAtSurfaceFromTo(ref int ANumber, double aStartX, double aEndX)
        {
            int i;
            double xcoor;
            double zcoor;
            for (i = FSurface.GetLowerBound(0) + 1; i <= FSurface.GetUpperBound(0); i++)
            {
                xcoor = FSurface[i].XCoor;
                zcoor = FSurface[i].ZCoor;
                if (xcoor > aStartX)
                {
                    if (xcoor < aEndX)
                    {
                        AddRestprofilePoint(ref ANumber, xcoor, zcoor);
                    }
                    else
                    {
                        break;
                    }
                }
            }
        }

        private bool CreateDOVRestprofile()
        {
            int LoopMax = 100;
            int LoopCount;
            double ZBot = DetermineLowestCoordinateFromLayers();
            int Lsign;
            double xLimit;
            int Number = 0;
            int LnDim;
            double LXMid;
            double LZMid;
            double xcoor;
            double zcoor;
            double ClaySlope;
            double SandSlope;
            bool LFound = false;
            bool result = false;
            bool LIsInAir = false;
            const double CCoorAccuracy = 0.001;
            ClaySlope = 2;
            SandSlope = 4;

            if ((FStabilitySide == TStabilitySideSet.LeftSide))
            {
                Lsign = -1;
            }
            else
            {
                Lsign = 1;
            }
            if ((Lsign > 0))
            {
                xLimit = DetermineMostRightCoordinateFromLayers();
            }
            else
            {
                xLimit = DetermineMostLeftCoordinateFromLayers();
            }

            //{ find possible itersection points with dike}
            List<double> lDikeIntersectionPoints = MStabDatafunctions.GetAllXAtSurface(FSurface, YStartDOV);
            LnDim = lDikeIntersectionPoints.Count;
            if (LnDim > 0)
            {
                if (FStabilitySide == TStabilitySideSet.RightSide)
                {
                    if (lDikeIntersectionPoints[0] < XStartDOV)
                    {
                        // there are a number of valid intersection points}
                        AddRestprofilePoint(ref Number, lDikeIntersectionPoints[0], YStartDOV);
                    }

                    for (int i = 1; i < lDikeIntersectionPoints.Count; i++)
                    {
                        if (lDikeIntersectionPoints[i] <= XStartDOV)
                        {
                            // check if line between two intersection points lies above surface
                            LXMid = 0.5*(lDikeIntersectionPoints[i] + lDikeIntersectionPoints[i - 1]);
                            LZMid = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, LXMid);
                            if (LZMid < YStartDOV)
                            {
                                // lies above surface so follow the surface}
                                addPointsAtSurfaceFromTo(ref Number, lDikeIntersectionPoints[i - 1],
                                                         lDikeIntersectionPoints[i]);
                            }
                            AddRestprofilePoint(ref Number, lDikeIntersectionPoints[i], YStartDOV);
                        }
                        else
                        {
                            if (lDikeIntersectionPoints[i - 1] < XStartDOV)
                            {
                                // former point before and this point after x meaning point lies in air or in dike
                                if (MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, XStartDOV) > YStartDOV)
                                {
                                    AddRestprofilePoint(ref Number, XStartDOV, YStartDOV);
                                    break;
                                }
                                else
                                {
                                    addPointsAtSurfaceFromTo(ref Number, lDikeIntersectionPoints[i - 1], YStartDOV);
                                    AddRestprofilePoint(ref Number, XStartDOV,
                                                        MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, XStartDOV));
                                    break;
                                }
                                ;
                            }
                        }
                    }
                }
            }

            if (Number == 0)
            {
                // geen punt gevonden hele talud langs lopen

                xcoor = FSurface[1].XCoor;
                zcoor = FSurface[1].ZCoor;

                AddRestprofilePoint(ref Number, xcoor, zcoor);
            }

            xcoor = FRestProfile[Number - 1].XCoor;
            zcoor = FRestProfile[Number - 1].YCoor;
            if (xcoor < XStartDOV)
            {
                addPointsAtSurfaceFromTo(ref Number, xcoor, XStartDOV);
            }
            double LTop = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, XStartDOV);
            if (YStartDOV < LTop)
            {
                AddRestprofilePoint(ref Number, XStartDOV, YStartDOV);
            }
            else
            {
                AddRestprofilePoint(ref Number, XStartDOV, LTop);
            }

            xcoor = FRestProfile[Number - 1].XCoor;
            zcoor = FRestProfile[Number - 1].YCoor;

            LoopCount = 0;
            do
            {
                LoopCount ++;
                FindIntersectionPoint(Lsign, ClaySlope, xLimit, ZBot, ref xcoor, ref zcoor, ref LFound, ref LIsInAir);
                if (LFound)
                {
                    // first see if the restprofile goes to clay

                    AddRestprofilePoint(ref Number, xcoor, zcoor);
                }
                if (!LFound && !LIsInAir)
                {
                    // If not so the restprofile passes sand
                    FindIntersectionPoint(Lsign, SandSlope, xLimit, ZBot, ref xcoor, ref zcoor, ref LFound, ref LIsInAir);
                    if (LFound)
                    {
                        AddRestprofilePoint(ref Number, xcoor, zcoor);
                    }
                }
                if (LIsInAir)
                {
                    // the dike is followed
                    MStabDatafunctions.GetnextpointOnSurface(FSurface, Lsign, ref xcoor, ref zcoor);
                    AddRestprofilePoint(ref Number, xcoor, zcoor);
                }
                if (!LFound && !LIsInAir)
                {
                    // Under a Clay slope sand is found while under a sand slope clay is found
                    // follow the boundary between those layers to the next point
                }
            } while (!((Math.Abs(xcoor - xLimit) < CCoorAccuracy) || (Math.Abs(zcoor - ZBot) < CCoorAccuracy) || (LoopCount >= LoopMax)));
            if ((LoopCount < LoopMax))
            {
                result = true;
            }
            // Bepaal ook de grensen van zone 3
            if ((Math.Abs(FXRoadLeft - FXRoadRight) > 0.1))
            {
                // er is een zone 3 aanwezig
                // zet de zone grensen eventuel goed
                if (FXRoadLeft > FXRoadRight)
                {
                    MStabDatafunctions.SwapDoubles(ref FXRoadLeft, ref FXRoadRight);
                }
                // Bereken de zone grensen
                FYRoadLeft = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, FXRoadLeft);
                FYRoadRight = MStabDatafunctions.GetZCoordinateSurfaceAtX(FSurface, FXRoadRight);
                FYBotLeft = ZBot;
                FYBotRight = ZBot;
                // Helling 1:1
                FXBotLeft = FXRoadLeft - Math.Abs(FYRoadLeft - FYBotLeft);
                FXBotRight = FXRoadRight - Math.Abs(FYRoadRight - FYBotLeft);
            }
            else
            {
                // geen zone 3, zet alle waarden op 0
                FYRoadLeft = 0;
                FYRoadRight = 0;
                FXBotLeft = 0;
                FXBotRight = 0;
                FYBotLeft = 0;
                FYBotRight = 0;
            }
            return result;
        }

        // =======================================================================================================================
        // Hier worden de punten van het restprofiel gemaakt.
        //
        // Date        ID   Modification
        // 2009-04-07       Created
        // =======================================================================================================================

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-23       Created
        // =======================================================================================================================
        private bool CirkelIntersectsRestProfile(double AXMid, double AZMid, double ARad, ref double AXRestProfile, ref double AZRestProfile)
        {
            bool result;
            int i = 0;
            int LIntersectionCount = 0;
            double xs1 = 0;
            double Ys1 = 0;
            double xs2 = 0;
            double Ys2 = 0;
            bool ready = false;

            i = 1;
            while ((i < FRestProfile.Count) && !ready)
            {
                MStabDatafunctions.Intersect_Circle_line(AXMid, AZMid, ARad, RestProfile[i - 1].XCoor, RestProfile[i].XCoor, RestProfile[i - 1].YCoor, RestProfile[i].YCoor, ref LIntersectionCount, ref xs1, ref Ys1, ref xs2, ref Ys2);
                if ((LIntersectionCount > 0))
                {
                    // in ieder geval een snijpunt kijk of het snijpunt op het lijnSTUK ligt
                    ready = Zones.PointBetweenPoints(RestProfile[i].XCoor, RestProfile[i].YCoor, RestProfile[i - 1].XCoor, RestProfile[i - 1].YCoor, xs1, Ys1);
                    if (ready)
                    {
                        AXRestProfile = xs1;
                        AZRestProfile = Ys1;
                    }
                    if ((LIntersectionCount > 1))
                    {
                        // twee snijpunten gevonden
                        if (Zones.PointBetweenPoints(RestProfile[i].XCoor, RestProfile[i].YCoor, RestProfile[i - 1].XCoor, RestProfile[i - 1].YCoor, xs2, Ys2))
                        {
                            // 2e snijpunt ligt op lijnSTUK
                            if (ready)
                            {
                                // beide punten op lijnstuk
                                if ((FStabilitySide == TStabilitySideSet.LeftSide))
                                {
                                    if ((xs2 > xs1))
                                    {
                                        AXRestProfile = xs2;
                                        AZRestProfile = Ys2;
                                    }
                                }
                                else
                                {
                                    if ((xs2 < xs1))
                                    {
                                        AXRestProfile = xs2;
                                        AZRestProfile = Ys2;
                                    }
                                }
                                // if ready
                            }
                            else
                            {
                                // alleen 2e punt ligt op lijnSTUK
                                ready = true;
                                AXRestProfile = xs2;
                                AZRestProfile = Ys2;
                            }
                        }
                    }
                }
                i ++;
            }
            result = ready;
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-07-08       Created
        // =======================================================================================================================
        private bool IsPointOnRestProfile(double AXCoor, double AZCoor)
        {
            bool result = false;
            int i = 0;
            bool Ready = false;
            double y3 = 0;
            double y4 = 0;
            double xs = 0;
            double ys = 0;
            i = 1;
            // kijken of intrede punt voor dijktafel ligt dan altijd zone 1
            if (FStabilitySide == TStabilitySideSet.LeftSide)
            {
                Ready = (AXCoor > RestProfile[0].XCoor);
            }
            else
            {
                Ready = (AXCoor < RestProfile[0].XCoor);
            }
            result = Ready;
            while ((i < FRestProfile.Count) && !Ready)
            {
                if ((AXCoor >= Math.Min(RestProfile[i].XCoor, RestProfile[i - 1].XCoor)) && (AXCoor <= Math.Max(RestProfile[i].XCoor, RestProfile[i - 1].XCoor)))
                {
                    Ready = true;
                    // maak en vertikale lijn door x
                    y3 = Math.Min(RestProfile[i].YCoor, RestProfile[i - 1].YCoor) - 10;
                    y4 = Math.Max(RestProfile[i].YCoor, RestProfile[i - 1].YCoor) + 10;
                    // Bepaal het snijpunt met het lijn
                    MStabDatafunctions.IntersectLines(RestProfile[i].XCoor, RestProfile[i].YCoor, RestProfile[i - 1].XCoor, RestProfile[i - 1].YCoor, AXCoor, y3, AXCoor, y4, ref xs, ref ys);
                    result = AZCoor < (ys + Constants.CGeoAccu);
                }
                i ++;
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-07-08       Created
        // =======================================================================================================================
        private bool CirkelIntersectsZone3(double Axmid, double AZMid, double ARad)
        {
            bool result = false;
            int LIntersectionCount = 0;
            double LXS1 = 0;
            double LXS2 = 0;
            double LZS1 = 0;
            double LZS2 = 0;
            LIntersectionCount = 0;
            if ((Math.Abs(FXRoadLeft - FXRoadRight) > Zones.CMinRoadWidth))
            {
                // see if circel intersects zone on left side
                MStabDatafunctions.Intersect_Circle_line(Axmid, AZMid, ARad, FXRoadLeft, FXBotLeft, FYRoadLeft, FYBotLeft, ref LIntersectionCount, ref LXS1, ref LZS1, ref LXS2, ref LZS2);
                if ((LIntersectionCount == 0))
                {
                    // If not may be circel intersects zone on right side
                    MStabDatafunctions.Intersect_Circle_line(Axmid, AZMid, ARad, FXRoadRight, FXBotRight, FYRoadRight, FYBotRight, ref LIntersectionCount, ref LXS1, ref LZS1, ref LXS2, ref LZS2);
                }
            }
            result = (LIntersectionCount > 0);
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2009-04-23       Created
        // =======================================================================================================================

        private bool DetermineXRestProfileAtY(double aYS, ref double aXS)
        {
            int i;
            double LX1;
            double LX2;
            double LXS = 0;
            double LYS = 0;
            bool result = false;

            LX1 = FRestProfile[0].XCoor;
            LX2 = FRestProfile[FRestProfile.Count - 1].XCoor;
            for (i = 1; i < FRestProfile.Count; i++)
            {
                if (MStabDatafunctions.IntersectLines(LX1, aYS, LX2, aYS, FRestProfile[i - 1].XCoor, FRestProfile[i - 1].YCoor,
                                                      FRestProfile[i].XCoor, FRestProfile[i].YCoor, ref LXS, ref LYS))
                {
                    result = true;
                    aXS = LXS;
                    break;
                }
            }
            return result;
        }

        private bool DetermineZSafeprofileAtX(double Ax, ref double Az)
        {
            int i;
            double XMax;
            double XMin;
            bool result = false;
            for (i = 0; i < SafeProfile.Count - 1; i++)
            {
                XMax = Math.Max(SafeProfile[i].XCoor, SafeProfile[i + 1].XCoor);
                XMin = Math.Min(SafeProfile[i].XCoor, SafeProfile[i + 1].XCoor);
                if ((Ax >= XMin) && (Ax <= XMax))
                {
                    Az = MStabDatafunctions.LinInpolY(SafeProfile[i].XCoor, SafeProfile[i].YCoor,
                                                      SafeProfile[i + 1].XCoor, SafeProfile[i + 1].YCoor, Ax);
                    result = true;
                    break;
                }
            }
            return result;
        }

        private bool LineIntersectsRestProfile(TPoints aCoorStart, TPoints aCoorEnd, ref double AXs, ref double AYs)
        {
            int i;
            bool LReady = false;
            //  LXS: Double;
            //  LYS: Double;

            i = 1;
            while ((i < FRestProfile.Count) && (! LReady))
            {
                LReady = MStabDatafunctions.IntersectLines(aCoorStart.XCoor, aCoorStart.ZCoor, aCoorEnd.XCoor, aCoorEnd.ZCoor,
                                                           FRestProfile[i - 1].XCoor, FRestProfile[i - 1].YCoor, FRestProfile[i].XCoor, FRestProfile[i].YCoor,
                                                           ref AXs, ref AYs);
                i++;
            }
            return LReady;
        }

        private bool PointInRestProfile(double XCoor, double YCoor)
        {
            double LXS = 0;
            bool result = false;

            if (FRestProfile.Count > 0)
            {
                if (DetermineXRestProfileAtY(YCoor, ref LXS))
                {
                    if ((FStabilitySide == TStabilitySideSet.LeftSide))
                    {
                        result = (XCoor > LXS);
                    }
                    else
                    {
                        result = (XCoor < LXS);
                    }
                }
            }
            return result;
        }
    }

    // end TZones
}

namespace Deltares.Stability.Calculation.Inner
{
    public class Zones
    {
        public const double CMinRoadWidth = 0.5;
        // =======================================================================================================================
        // Point is calculated on intersection
        // Check if it's indeed on the actual line
        //
        // Date        ID   Modification
        // 2009-04-23       Created
        // =======================================================================================================================
        public static bool PointBetweenPoints(double AX1, double AY1, double AX2, double AY2, double AXS, double AYS)
        {
            bool result;
            result = (AXS >= Math.Min(AX1, AX2)) && (AXS <= Math.Max(AX1, AX2)) && (AYS >= Math.Min(AY1, AY2)) && (AYS <= Math.Max(AY1, AY2));
            return result;
        }
    }

    // end Zones
}