// Class name: TNail
//
// Description:
//
// Copyright (c) 1999-2012 Deltares
//
// Date        ID      Modification
// 2012-07-24  Berg_pl Created
// ===========================================================================================

using System;
using System.Collections.Generic;
using Deltares.Mathematics;
using Deltares.Standard;

namespace Deltares.Stability.Calculation.Inner
{
    public class TNail
    {
        #region Variable: Constants

        private enum TIntersectType
        {
            TiIntersect,
            TiCrit3,
            TiCrit4,
            TiUnexpected
        };

        private const double CCoorAccuracy = 0.001;
        private const double D2R = Math.PI/180;

        #endregion

        #region Variable: Calculated values

        private double FAnchorLength;
        private double FCrit1;
        private double FCrit2;
        private double FDiameter;
        private double FFDMax;
        private double FFNMax;

        //FPointGeometry: TDGeometry;
        //FPointSoilList: TDSoilList;

        private double FL0;
        private int FLayerIndex;
        private double FPU;
        private double FResistingMoment;
        private List<Point2D> FSafetyArea = new List<Point2D>();
        private TStabSoilData FSoil;
        private double FXEnd;
        private double FXIntersection;
        private double FYEnd;
        private double FYIntersection;

        #endregion

        #region Variable: Properties

        public double ResistingMoment
        {
            get
            {
                return FResistingMoment;
            }
            set
            {
                FResistingMoment = value;
            }
        }

        public string NailName
        {
            get
            {
                return FNailName;
            }
            set
            {
                FNailName = value;
            }
        }

        public double X
        {
            get
            {
                return FXStart;
            }
            set
            {
                FXStart = value;
            }
        }

        public double Y
        {
            get
            {
                return FYStart;
            }
            set
            {
                FYStart = value;
            }
        }

        public double DirectionAngle
        {
            get
            {
                return FAngle;
            }
            set
            {
                FAngle = value;
            }
        }

        public double CriticalAngle
        {
            get
            {
                return FCriticalAngle;
            }
            set
            {
                FCriticalAngle = value;
            }
        }

        public double HorizontalSpacing
        {
            get
            {
                return FHorizontalSpacing;
            }
            set
            {
                FHorizontalSpacing = value;
            }
        }

        public double Length
        {
            get
            {
                return FLength;
            }
            set
            {
                FLength = value;
            }
        }

        public double DiameterNail
        {
            get
            {
                return FDiameterNail;
            }
            set
            {
                FDiameterNail = value;
            }
        }

        public double DiameterGrout
        {
            get
            {
                return FDiameterGrout;
            }
            set
            {
                FDiameterGrout = value;
            }
        }

        public double YieldForce
        {
            get
            {
                return FYieldForce;
            }
            set
            {
                FYieldForce = value;
            }
        }

        public double PlasticMoment
        {
            get
            {
                return FPlasticMoment;
            }
            set
            {
                FPlasticMoment = value;
            }
        }

        public double BendingStiffness
        {
            get
            {
                return FBendingStiffness;
            }
            set
            {
                FBendingStiffness = value;
            }
        }

        public bool IsFacingUsed
        {
            get
            {
                return FIsFacingUsed;
            }
            set
            {
                FIsFacingUsed = value;
            }
        }

        public TTauMAxCalcTypeSet TaumaxType
        {
            get
            {
                return FTaumaxType;
            }
            set
            {
                FTaumaxType = value;
            }
        }

        public bool UseShearstressTable
        {
            get
            {
                return FUseShearstressTable;
            }
            set
            {
                FUseShearstressTable = value;
            }
        }

        public bool UseLateralStressTable
        {
            get
            {
                return FUseLateralStressTable;
            }
            set
            {
                FUseLateralStressTable = value;
            }
        }

        public double MaxPullForce
        {
            get
            {
                return FMaxPullForce;
            }
            set
            {
                FMaxPullForce = value;
            }
        }

        public double MaxMom
        {
            get
            {
                return FMaxMom;
            }
            set
            {
                FMaxMom = value;
            }
        }

        public double EI
        {
            get
            {
                return FEI;
            }
            set
            {
                FEI = value;
            }
        }

        public TShearStressTable ShearStressTable
        {
            get
            {
                return FShearStressTable;
            }
            set
            {
                FShearStressTable = value;
            }
        }

        public TLateralStressTable LateralStressCurves
        {
            get
            {
                return FLateralStressCurves;
            }
            set
            {
                FLateralStressCurves = value;
            }
        }

        #endregion

        private double FAngle;
        private double FBendingStiffness;
        private double FCriticalAngle;
        private double FDiameterGrout;
        private double FDiameterNail;
        private double FEI;
        private TPoints[] FFreaticLine;
        private double FGammaWater;
        private double FHorizontalSpacing;
        private TInterfaceCalls FInterfaceCalls;
        private TInterfaceData FInterfaceData;
        private bool FIsFacingUsed;
        private TLateralStressTable FLateralStressCurves;
        private double FLength;

        private double FMaxMom;
        private double FMaxPullForce;
        private string FNailName = String.Empty;
        private TPoints[][] FPLLines;
        private double FPlasticMoment;

        private TShearStressTable FShearStressTable;
        private TStabSoilData[] FSoilData;

        private TPoints[] FSurfaceLine;
        private TTauMAxCalcTypeSet FTaumaxType;
        private bool FUseLateralStressTable;
        private bool FUseShearstressTable;
        private TWaterLineRecord[] FWaterLines;
        private double FXStart;
        private double FYStart;
        private double FYieldForce;

        public TNail(TInterfaceData aStabilityObjects)
        {
            FInterfaceData = aStabilityObjects;
            InitNail();
        }

        public TNail(TNail aNail)
        {
            InitNail();

            if (aNail != null)
            {
                FXStart = aNail.X;
                FYStart = aNail.Y;
                FLength = aNail.Length;
                FAngle = aNail.DirectionAngle;
                FIsFacingUsed = aNail.IsFacingUsed;
                //new
                FTaumaxType = aNail.TaumaxType;
                FCriticalAngle = aNail.CriticalAngle;
                FDiameterNail = aNail.DiameterNail;
                FDiameterGrout = aNail.DiameterGrout;
                FEI = aNail.EI;
                FHorizontalSpacing = aNail.HorizontalSpacing;
                FMaxPullForce = aNail.MaxPullForce;
                FMaxMom = aNail.MaxMom;
                FLateralStressCurves.Copy(aNail.LateralStressCurves);
                FShearStressTable.Copy(aNail.ShearStressTable);
                FUseShearstressTable = aNail.UseShearstressTable;
                FUseLateralStressTable = aNail.UseLateralStressTable;

                FResistingMoment = aNail.ResistingMoment;
            }
        }

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

        public double CalculateResistingMoment(double AXMid, double AYMid, double ARad, TInterfaceData aInterfaceData)
        {
            double LFNail;
            double LCentreAngle;
            double LParallelAngle;
            double LForceParallel;
            double LForcePerpen;
            double LCrit1;
            double LMomFN;
            double LMomFD;
            double LForceAngle;
            double LNeedleAngle;
            double LDifAngle;
            double LInitAngle;
            FInterfaceData = aInterfaceData;
            FInterfaceCalls = aInterfaceData.InterfaceCalls;
            FFNMax = 0;
            FFDMax = 0;

            FResistingMoment = 0;

            try
            {
                if (FindIntersectionPoints(AXMid, AYMid, ARad))
                {
                    if (CalculateNeededParameters())
                    {
                        FCrit1 = CalculateFNMax();

                        FCrit2 = 0.5*FPU*FDiameter*FL0;

                        LCentreAngle = DArctan2(AYMid - FYIntersection, AXMid - FXIntersection);
                        LNeedleAngle = FAngle*Math.PI/180 + Math.PI;
                        if (Math.Abs(LNeedleAngle - 2*Math.PI) < 0.00001)
                        {
                            LNeedleAngle = 0;
                        }

                        LDifAngle = (LNeedleAngle - LCentreAngle);
                        LInitAngle = Math.PI/2 - Math.Abs(LDifAngle);

                        if (Math.Abs(LInitAngle) <= FCriticalAngle)
                        {
                            // Only FN 
                            LCrit1 = CalculateFNMax();
                            FFNMax = Math.Min(LCrit1, FMaxPullForce);
                            FFDMax = 0;
                        }

                        else
                        {
                            if (Math.Abs(LInitAngle) >= Math.PI/2 - FCriticalAngle)
                            {
                                FFNMax = 0;
                                // only FD 
                                FFDMax = CalculateFdMaxShearOnly();
                            }
                            else
                            {
                                CalculateFdAndFn(LInitAngle);
                            }
                        }
                        LFNail = Math.Sqrt(FFNMax*FFNMax + FFDMax*FFDMax);
                        LForceAngle = DArctan(FFDMax, FFNMax);
                        LParallelAngle = LInitAngle - LForceAngle;

                        LForceParallel = LFNail*Math.Cos(LParallelAngle);
                        LForcePerpen = LFNail*Math.Sin(LParallelAngle);

                        // Momenten bepalen 
                        LMomFN = Math.Abs(LForcePerpen)*Math.Tan(FSoil.Phi*Math.PI/180)*ARad;
                        LMomFD = Math.Abs(LForceParallel)*ARad;

                        FResistingMoment = LMomFN + LMomFD;
                        if (FHorizontalSpacing > 0)
                        {
                            FResistingMoment = FResistingMoment/FHorizontalSpacing;
                        }
                    }
                }
            }
            finally {}
            return FResistingMoment;
        }

        private void InitNail()
        {
            FNailName = "NoName";
            FXStart = 0.0;
            FYStart = 0.0;
            FAngle = 0.0;
            FHorizontalSpacing = 0;
            FLength = 0;
            FDiameterNail = 0;
            FDiameterGrout = 0;
            FPlasticMoment = 0;
            FBendingStiffness = 0;
            //TODO: set to tabel for now, fix import; does not work yet
            FTaumaxType = TTauMAxCalcTypeSet.tmTabel;
            FIsFacingUsed = false;
            FUseShearstressTable = false;
            FUseLateralStressTable = false;
            FMaxPullForce = 0;
            FMaxMom = 0;
            FInterfaceCalls = new TInterfaceCalls();
            FSoilData = FInterfaceData.StabSoilData;
            if (FInterfaceData.CurrentWaternet.PhreaticLineNumber >= 0)
            {
                FFreaticLine = FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber];
            }
            FSurfaceLine = FInterfaceData.SurfaceLine;
            FPLLines = FInterfaceData.CurrentWaternet.PLLines;
            FWaterLines = FInterfaceData.CurrentWaternet.Waterlines;
            FGammaWater = FInterfaceData.CurrentWaternet.GammaWater;
        }

        private bool FindIntersectionPoints(double AXMid, double AYMid, double ARad)
        {
            bool result = false;

            //Find the endpoints of the nail
            FXEnd = FXStart + FLength*Math.Cos(FAngle*D2R);
            FYEnd = FYStart + FLength*Math.Sin(FAngle*D2R);

            //find intersection point(s) with circle
            var intersections = Routines2D.IntersectCircleline(AXMid, AYMid, ARad,
                                                               FXStart, FXEnd, FYStart, FYEnd);
            if (intersections.Count == 1)
            {
                FXIntersection = intersections[0].X;
                FYIntersection = intersections[0].Y;
                result = true;
            }
            else if (intersections.Count > 1)
            {
                FXIntersection = intersections[1].X;
                FYIntersection = intersections[1].Y;
                result = true;
            }
            return result;
        }

        private TStabSoilData DetermineSoilFromCoordinates(double Ax, double Ay)
        {
            int FLayerIndex = FInterfaceCalls.DetermineLayerNumberFromPoint(Ax, Ay);
            return FSoilData[FLayerIndex];
        }

        private double DetemineEffectiveStress(double Ax, double Ay)
        {
            var LVertical = new TVertical();
            LVertical.SoilData.AddRange(FSoilData);
            LVertical.XCoordinate = Ax;
            LVertical.ZCoordinate = Ay;
            LVertical.PLLines = FPLLines;
            LVertical.GammaWater = FGammaWater;
            LVertical.WaterLines.AddRange(FWaterLines);
            LVertical.FreaticLine = FFreaticLine;
            LVertical.SurfaceLine = FSurfaceLine;
            LVertical.FillVerticalWithWaterData();
            LVertical.DegreeOfConsolidationSoil = FInterfaceData.DegreeOfConsolidationSoil;
            LVertical.InterfaceCalls = FInterfaceCalls;
            LVertical.FillVertical();
            return LVertical.DetermineEffectiveStress(Ay);
        }

        private double DeteminePreloadEffectiveStress(double Ax, double AYBot, double AYTop)
        {
            var LVertical = new TVertical();
            LVertical.SoilData.AddRange(FSoilData);
            LVertical.XCoordinate = Ax;
            LVertical.ZCoordinate = AYBot;
            LVertical.PLLines = FPLLines;
            LVertical.GammaWater = FGammaWater;
            LVertical.WaterLines.AddRange(FWaterLines);
            LVertical.FreaticLine = FFreaticLine;
            LVertical.SurfaceLine = FSurfaceLine;
            LVertical.FillVerticalWithWaterData();
            LVertical.DegreeOfConsolidationSoil = FInterfaceData.DegreeOfConsolidationSoil;
            LVertical.InterfaceCalls = FInterfaceCalls;
            LVertical.FillVertical();
            return LVertical.DeterminePreLoadEffectiveStress(AYBot, AYTop);
        }

        private bool CalculateNeededParameters()
        {
            bool result = false;
            double LEs;
            double LStress;

            FLayerIndex = FInterfaceCalls.DetermineLayerNumberFromPoint(FXIntersection, FYIntersection);
            FDiameter = Math.Max(FDiameterNail, FDiameterGrout);

            {
                FSoil = FSoilData[FLayerIndex];
                LStress = DetemineEffectiveStress(FXIntersection, FYIntersection);
                FPU = GetPu(FXIntersection, FYIntersection, LStress);
                LEs = GetEs(LStress);
                LEs = Math.Max(LEs, 1);

                FL0 = Math.Pow(4*FEI/LEs, 0.25);
                if (FL0 <= 0.00001)
                {
                    FL0 = 0.00001;
                }
                result = true;
            }

            return result;
        }

        private double GetPu(double AXCoor, double AYCoor, double AStress)
        {
            double LKp;
            double LSin;
            double LSchelp;
            double LCu;
            double LPu;
            double LLength;
            TMatStrengthTypeSet LStrength;

            LPu = 0;

            if (FUseLateralStressTable)
            {
                LLength = Routines2D.Compute2DDistance(AXCoor, AYCoor, FXStart, FYStart);
                LPu = FLateralStressCurves.RetrieveYatX(LLength);
            }
            else
            {
                LStrength = FSoil.MaterialStrength;
                if (LStrength == TMatStrengthTypeSet.mstNone)
                {
                    LStrength = FInterfaceData.ModelData.DefaultStrength;
                }

                if ((LStrength == TMatStrengthTypeSet.mstCphi) ||
                    (LStrength == TMatStrengthTypeSet.mstStressTab) ||
                    (LStrength == TMatStrengthTypeSet.mstPseudoStressTab))
                {
                    // LSchelp = LSoil.SchelpFactor;
                    // ongedraineerd aangenomen 
                    LSchelp = 3;
                    // non cohesive material
                    LSin = Math.Sin(FSoil.Phi*Math.PI/180);
                    LKp = (1 + LSin)/(1 - LSin);
                    LPu = LSchelp*AStress*LKp;
                }

                if ((LStrength == TMatStrengthTypeSet.mstMeasuredCu)
                    || (LStrength == TMatStrengthTypeSet.mstCalculatedCu)
                    || (LStrength == TMatStrengthTypeSet.mstGradientCu))
                {
                    LCu = GetCuValuesAtCoor(FXIntersection, FYIntersection);
                    LPu = 9*LCu;
                }
            }

            return LPu;
        }

        private double GetEs(double AStress)
        {
            double LAlfa;

            double result = 1;
            LAlfa = GetAlfaFromSoilType(FSoil);
            result = (AStress*LAlfa)/(FSoil.RRatio*((2.0/9.0)*Math.Pow(2.65, LAlfa) + (LAlfa/6.0)));
            return result;
        }

        private double GetAlfaFromSoilType(TStabSoilData aSoil)
        {
            double LAlfa = 1;
            switch (aSoil.SoilType)
            {
                case TSoilTypeSet.stSand:
                    LAlfa = 1.0/3.0;
                    break;
                case TSoilTypeSet.stLoam:
                    LAlfa = 1.0/2.0;
                    break;
                case TSoilTypeSet.stClay:
                    LAlfa = 2.0/3.0;
                    break;
                case TSoilTypeSet.stPeat:
                    LAlfa = 1.0;
                    break;
                case TSoilTypeSet.stGravel:
                    LAlfa = 1/4;
                    break;
            }
            return LAlfa;
        }

        private double CalculateFNMax()
        {
            double LTauMax = 0;
            double LXM;
            double LYM;

            if (FIsFacingUsed)
            {
                FAnchorLength = Routines2D.Compute2DDistance(FXIntersection, FYIntersection, FXEnd, FYEnd);
                LXM = 0.5*(FXIntersection + FXEnd);
                LYM = 0.5*(FYIntersection + FYEnd);
                LTauMax = GetTauMax(LXM, LYM);
            }
            else
            {
                double L1 = Routines2D.Compute2DDistance(FXIntersection, FYIntersection, FXEnd, FYEnd);
                double L2 = Routines2D.Compute2DDistance(FXIntersection, FYIntersection, FXStart, FYStart);
                if (L1 < L2)
                {
                    FAnchorLength = L1;
                    LXM = 0.5*(FXIntersection + FXEnd);
                    LYM = 0.5*(FYIntersection + FYEnd);
                    LTauMax = GetTauMax(LXM, LYM);
                }
                else
                {
                    FAnchorLength = L2;
                    LXM = 0.5*(FXIntersection + FXStart);
                    LYM = 0.5*(FYIntersection + FYStart);
                    LTauMax = GetTauMax(LXM, LYM);
                }
            }

            return LTauMax*FDiameter*FAnchorLength*Math.PI;
        }

        private void CalculateFdAndFn(double AInitAngle)
        {
            int i;
            var LXArr = new List<double>();
            double LXE1;
            double LXE3;
            double LXE4;
            double LXLast;
            double LX;
            double LY = 0.0;
            double LAngle;
            double LTan;
            double LDist;
            double LMax;
            int LMaxNum;
            TIntersectType LIntersection;

            // calculate criteria values on FN axes
            LXE1 = FCrit1;
            LXE3 = Criterium3Invert(0);
            LXE4 = Criterium4Invert(0);
            LXLast = LXE1;

            if (LXE3 < LXLast)
            {
                LXLast = LXE3;
            }

            if (LXE4 < LXLast)
            {
                LXLast = LXE4;
            }

            LXArr.Insert(0, 0.0);
            LXArr.Insert(1, LXLast);

            // voeg eventueel de knopen toe
            AddXCoor(LXArr, LXE1);
            AddXCoor(LXArr, LXE3);
            AddXCoor(LXArr, LXE4);

            // kijk of er meer knopen zijn
            // snijpunt Criterium 2 em 3
            LXE3 = Criterium3(FCrit2);
            AddXCoor(LXArr, LXE3);

            // snijpunt Criterium 2 em 4
            LXE4 = Criterium4(FCrit2);
            AddXCoor(LXArr, LXE4);

            // snijpunt Criterium 3 em 4
            var result = CalcIntersectionCrit3AndCrit4();
            LXE4 = result.Item2; //LXE4
            LY = result.Item3; // LY
            LIntersection = result.Item1; //intersection type
            if (LIntersection == TIntersectType.TiIntersect)
            {
                AddXCoor(LXArr, LXE4);
            }

            // Nu zijn alle x coordinaten (Fn) waarop wat gebeurt bekend Zoek nu per punt de laagste bijbehorende FD
            for (i = 0; i < LXArr.Count; i++)
            {
                LX = LXArr[i];
                LY = GetLowestFd(LX);
                FSafetyArea.Add(new Point2D(LX, LY));

                // if point is criterium 1 then ad vertical point
                if (Math.Abs(FCrit1 - LX) < CCoorAccuracy)
                {
                    FSafetyArea.Add(new Point2D(LX, 0));
                }
            }

            // nu het veiligheidsdomein is gevonden de grootste kracht evenwijdig aan de cirkel opzoeken
            LAngle = Math.PI/2 - Math.Abs(AInitAngle);
            LTan = Math.Tan(LAngle);
            LMax = 0;
            LMaxNum = 0;
            for (i = 0; i < FSafetyArea.Count; i++)
            {
                LDist = FSafetyArea[i].X + FSafetyArea[i].Y/LTan;
                if (LDist > LMax)
                {
                    LMax = LDist;
                    LMaxNum = i;
                }
            }

            FFNMax = FSafetyArea[LMaxNum].X;
            FFDMax = FSafetyArea[LMaxNum].Y;
        }

        private void AddXCoor(List<double> AXarr, double AX)
        {
            int LPos;
            int i;
            bool LFound = false;

            if ((AX < AXarr[AXarr.Count - 1] - CCoorAccuracy) && (AX > 0))
            {
                // punt valt binnen bereik
                for (i = 1; i <= AXarr.Count; i++)
                {
                    if (Math.Abs(AX - AXarr[i - 1]) < CCoorAccuracy)
                    {
                        LFound = true;
                    }
                }

                if (!LFound)
                {
                    // punt tussenvoegen
                    i = 1;
                    while ((AX < AXarr[i - 1]) && (i < AXarr.Count))
                    {
                        i++;
                    }

                    LPos = i;

                    //ANDim++;
                    //AXarr.Capacity = ANDim;
                    for (i = AXarr.Count - 1; i > LPos; i--)
                    {
                        AXarr.Insert(i, AXarr[i - 1]);
                        //AXarr[i] = AXarr[i - 1];
                    }
                    AXarr.Insert(LPos, AX);
                    //AXarr[LPos] = AX;
                }
            }
        }

        private double GetLowestFd(double ax)
        {
            double result = FCrit2;

            double lfd = Criterium3(ax);
            if (lfd < result)
            {
                result = lfd;
            }

            lfd = Criterium4(ax);
            if (lfd < result)
            {
                result = lfd;
            }

            return result;
        }

        private double Criterium3(double afn)
        {
            double result = 0;
            if (FMaxPullForce > afn)
            {
                result = 0.5*Math.Sqrt(FMaxPullForce*FMaxPullForce - afn*afn);
            }
            return result;
        }

        private double Criterium3Invert(double AFD)
        {
            double result = 0;
            if (FMaxPullForce > 2*AFD)
            {
                result = Math.Sqrt(FMaxPullForce*FMaxPullForce - 4*AFD*AFD);
            }
            return result;
        }

        private double Criterium4(double afn)
        {
            double result = FMaxPullForce;

            if (FMaxPullForce > 0)
            {
                double LRatio = afn/FMaxPullForce;
                if (FMaxPullForce > afn)
                {
                    result = 1.62*FMaxMom/FL0*(1 - LRatio*LRatio) + 0.24*FPU*FDiameter*FL0;
                }
            }
            return result;
        }

        private double Criterium4Invert(double AFD)
        {
            double result = 0;

            double LFactor1 = 1.62*FMaxMom/FL0;
            double LFactor2 = 0.24*FPU*FDiameter*FL0;
            if (AFD < LFactor1 + LFactor2)
            {
                double LRatio = FMaxPullForce*FMaxPullForce/LFactor1;
                result = -LRatio*(AFD - LFactor1 - LFactor2);
                result = Math.Sqrt(result);
            }
            return result;
        }

        private Tuple<TIntersectType, double, double> CalcIntersectionCrit3AndCrit4()
        {
            double LFactor1;
            double LFactor2;
            double La;
            double Lb;
            double Lc;
            double LDiscriminant;
            double Ld3;
            double Ld4;
            double Ln3;
            double Ln4;
            double LY = 0;
            double LXE4 = 0;
            var intersectType = TIntersectType.TiUnexpected;
            try
            {
                LFactor1 = 1.62*FMaxMom/FL0;
                LFactor2 = 0.24*FPU*FDiameter*FL0;

                // both criteria together in a quadratic equation of AFD (a AFD2 + b * AFd + c = 0
                La = 4*LFactor1/(FMaxPullForce*FMaxPullForce);
                Lb = -1;
                Lc = LFactor2;
                LDiscriminant = Lb*Lb - 4*La*Lc;
                if (LDiscriminant > 0)
                {
                    LY = (-Lb - Math.Sqrt(LDiscriminant))/(2*La);
                    LXE4 = Math.Sqrt(FMaxPullForce*FMaxPullForce - 4*LY*LY);
                    intersectType = TIntersectType.TiIntersect;
                }
                else
                {
                    // kijk of het ene criterium binnen het andere ligt
                    Ln3 = Criterium3Invert(0);
                    Ln4 = Criterium4Invert(0);
                    Ld3 = Criterium3(0);
                    Ld4 = Criterium4(0);
                    if ((Ln3 <= Ln4) && (Ld3 <= Ld4))
                    {
                        intersectType = TIntersectType.TiCrit3;
                    }
                    if ((Ln4 <= Ln3) && (Ld4 <= Ld3))
                    {
                        intersectType = TIntersectType.TiCrit4;
                    }
                }
            }
            catch
            {
                LY = 0;
                LXE4 = 0;
                intersectType = TIntersectType.TiUnexpected;
            }
            var result = new Tuple<TIntersectType, double, double>(intersectType, LXE4, LY);
            return result;
        }

        private double CalculateFdMaxShearOnly()
        {
            //  LCu: Double;
            // TODO - oHan : Wat als met methode c, phi gerekend wordt.Moet de Cu uit het materiaal gehaald worden of uit cu berekend, gemeten if ratio? 

            double LCrit2 = 0.5*FPU*FDiameter*FL0;
            double LCrit4 = 1.62*FMaxMom/FL0 + 0.24*FPU*FDiameter*FL0;

            double result = Math.Min(FMaxPullForce/2, LCrit2);
            result = Math.Min(result, LCrit4);
            return result;
        }

        private double GetTauMax(double Ax, double Ay)
        {
            int LLayerIndex;
            double LYSurface;
            double LDepth;
            TStabSoilData LSoil; // pointer
            double LFactor;
            double LAlfaT;
            double Lqc;
            double LStress;
            double LLength;

            double result = 0;

            // tijdelijk voor testen
            LFactor = 1;
            LLayerIndex = FInterfaceCalls.DetermineLayerNumberFromPoint(Ax, Ay);
            if (LLayerIndex >= 0)
            {
                LSoil = FSoilData[FLayerIndex];

                if (FUseShearstressTable)
                {
                    LLength = Routines2D.Compute2DDistance(Ax, Ay, FXStart, FYStart);
                    result = FShearStressTable.RetrieveYatX(LLength);
                }
                else
                {
                    LYSurface = MStabDatafunctions.ZTopAtSurface(Ax, FSurfaceLine);
                    LDepth = LYSurface - Ay;

                    switch (FTaumaxType)
                    {
                            // Via tabel
                        case TTauMAxCalcTypeSet.tmTabel:
                            {
                                LStress = DetemineEffectiveStress(Ax, Ay);
                                result = LSoil.BondStressTable.GetTauFromStress(LStress);
                                break;
                            }
                        case TTauMAxCalcTypeSet.tmGiven:
                            result = LFactor*LSoil.CuTop;
                            break;
                        case TTauMAxCalcTypeSet.tmCPT:
                            {
                                LAlfaT = 0.01;
                                //  via cpt values
                                Lqc = GetQcAtDepth(LDepth);
                                if (LSoil.SoilType == TSoilTypeSet.stSand)
                                {
                                    LAlfaT = 0.012;
                                }

                                if (LSoil.SoilType == TSoilTypeSet.stClay)
                                {
                                    LAlfaT = 0.025;
                                    if (Lqc < 1)
                                    {
                                        if (LDepth/FDiameter < 20)
                                        {
                                            LAlfaT = 0.02;
                                        }
                                    }
                                }
                                result = LAlfaT*Lqc;
                                break;
                            }
                    }
                }
            }
            return result;
        }

        private double GetQcAtDepth(double aDepth)
        {
            if (aDepth < 10)
            {
                return 1;
            }

            // TODO : replace with a real function
            return aDepth < 20 ? 2 : 3;
        }

        private double GetCuValuesAtCoor(double aXCoor, double aYCoor)
        {
            double ygb;
            double ygo;
            double LInterpolationFactorCu;
            double LPreLoadStress;
            double LUsedStress;
            double LPop;
            double LEffStress;
            TStabSoilData LSoil;

            double result = 0;
            LSoil = DetermineSoilFromCoordinates(aXCoor, aYCoor);
            //do begin LPop = LSoil .PopValueUsed;
            LPop = LSoil.POP;
            ygb = FInterfaceCalls.GetZTopFromLayer(aXCoor, aYCoor);
            ygo = FInterfaceCalls.GetZBottomFromLayer(aXCoor, aYCoor);

            // ... Cu Measured ...
            if (LSoil.MaterialStrength == TMatStrengthTypeSet.mstMeasuredCu)
            {
                if (ygb > ygo)
                {
                    LInterpolationFactorCu = (aYCoor - ygo)/(ygb - ygo);
                }
                else
                {
                    LInterpolationFactorCu = 1;
                }

                return LInterpolationFactorCu*(LSoil.CuTop - LSoil.CuBottom) + LSoil.CuBottom;
            }

            // ... Cu Gradient (not yet probabilistic) ...
            if (LSoil.MaterialStrength == TMatStrengthTypeSet.mstGradientCu)
            {
                return (LSoil.CuTop + (ygb - aYCoor)*LSoil.CuGradient);
            }

            // ... Cu calculated ...
            if (LSoil.MaterialStrength == TMatStrengthTypeSet.mstCalculatedCu)
            {
                LEffStress = DetemineEffectiveStress(aXCoor, aYCoor);
                double PreloadZ = FInterfaceCalls.DetermineZReferencelevelCuCalculated(aXCoor);
                double surfaceZ = MStabDatafunctions.ZTopAtSurface(aXCoor, FSurfaceLine);
                if (PreloadZ != surfaceZ)
                {
                    LPreLoadStress = DeteminePreloadEffectiveStress(aXCoor, aYCoor, PreloadZ);
                }
                else
                {
                    LPreLoadStress = LEffStress;
                }

                LUsedStress = Math.Max(LPreLoadStress + LPop, LEffStress);
                return ((LUsedStress)*LSoil.RatioCuPc);
            }
            return result;
        }

        private double DArctan(double y, double x)
        {
            double result;
            if (Math.Abs(x - 0) < CCoorAccuracy)
            {
                if (y < 0)
                {
                    result = -Math.PI/2;
                }
                else
                {
                    result = Math.PI/2;
                }
            }
            else
            {
                result = Math.Atan(y/x);
            }
            return result;
        }

        private double DArctan2(double y, double x)
        {
            double result;
            if (Math.Abs(x - 0) < CCoorAccuracy)
            {
                if (y < 0)
                {
                    result = -Math.PI/2;
                }
                else
                {
                    result = Math.PI/2;
                }
            }
            else
            {
                result = Math.Atan2(y, x);
            }
            return result;
        }
    }
}