using System;
using System.Collections;
using System.Diagnostics;
using Deltares.Mathematics;

// =======================================================================================================================
// Class name:
//
// Description:
//
// Copyright (c) 2001-2010 Deltares
//
// Date        ID   Modification
// 2008-04-23  Best Created
//
// (Probabilistic voorlopig niet)
// =======================================================================================================================

namespace Deltares.Stability.Calculation.Inner
{
    public class TSliceFill
    {
        private bool FAreNegativeValuesAllowed = false;
        private TProbInputValuesSet FDefaultInputValues;
        private TMatStrengthTypeSet FDefaultStrength;
        private double[] FDegreeOfConsolidationEarthQuake;
        private double[] FDegreeOfConsolidationLoads;
        private double[][] FDegreeOfConsolidationSoil;
        private TEarthquakeData FEarthquakeData;
        private TPoints[] FFreaticLine;
        private double FGammaWater = 0;
        private TInterfaceCalls FInterfaceCalls = null;
        private bool FIsAlternativeStrengthUsed = false;
        private bool FIsCapillairLoadUsed = false;
        private bool FIsGDInternalVersion = false;
        private bool FIsGlobalMesurement = false;
        private bool FIsMSeepnetUsed = false;
        private bool FIsProbabilisticCalculation = false;
        private bool FIsUpliftCalculation = false;
        private int FLeftNumber = 0;
        private TLineLoad[] FLineLoad;
        private string FMSeepfilename = String.Empty;
        private int FNSoils = 0;
        private TPoints[][] FPLLines;
        private TPoints[] FPreconStressLine;
        private TProbCalculationType FProbabilisticCalculationType;
        private int FRightNumber = 0;
        private TStabSoilData[] FSoilData;
        private TStabilitySideSet FStabilitySide;
        private TStressTable[] FStressTable;
        private TPoints[] FSurfaceLine;
        private TUniformLoad[] FUniformLoad;
        private TWaterLineRecord[] FWaterLines;

        public TSliceFill()
            : base()
        {
            // FInterfaceCalls := TInterfaceCalls.Create;
        }

        public TMatStrengthTypeSet DefaultStrength
        {
            get
            {
                return FDefaultStrength;
            }
            set
            {
                FDefaultStrength = value;
            }
        }

        public TProbInputValuesSet DefaultInputValues
        {
            get
            {
                return FDefaultInputValues;
            }
            set
            {
                FDefaultInputValues = value;
            }
        }

        public TProbCalculationType ProbabilisticCalculationType
        {
            get
            {
                return FProbabilisticCalculationType;
            }
            set
            {
                FProbabilisticCalculationType = value;
            }
        }

        public TStabSoilData[] SoilData
        {
            get
            {
                return FSoilData;
            }
            set
            {
                FSoilData = value;
            }
        }

        public TPoints[] FreaticLine
        {
            get
            {
                return FFreaticLine;
            }
            set
            {
                FFreaticLine = value;
            }
        }

        public TPoints[] PreconStressLine
        {
            get
            {
                return FPreconStressLine;
            }
            set
            {
                FPreconStressLine = value;
            }
        }

        public TPoints[] SurfaceLine
        {
            get
            {
                return FSurfaceLine;
            }
            set
            {
                FSurfaceLine = value;
            }
        }

        public bool IsMSeepNetUsed
        {
            get
            {
                return FIsMSeepnetUsed;
            }
            set
            {
                FIsMSeepnetUsed = value;
            }
        }

        public bool IsCapillairLoadUsed
        {
            get
            {
                return FIsCapillairLoadUsed;
            }
            set
            {
                FIsCapillairLoadUsed = value;
            }
        }

        public bool AreNegativeValuesAllowed
        {
            get
            {
                return FAreNegativeValuesAllowed;
            }
            set
            {
                FAreNegativeValuesAllowed = value;
            }
        }

        public string MSeepfilename
        {
            get
            {
                return FMSeepfilename;
            }
            set
            {
                FMSeepfilename = value;
            }
        }

        public bool IsUpliftCalculation
        {
            get
            {
                return FIsUpliftCalculation;
            }
            set
            {
                FIsUpliftCalculation = value;
            }
        }

        public bool IsGlobalMesurement
        {
            get
            {
                return FIsGlobalMesurement;
            }
            set
            {
                FIsGlobalMesurement = value;
            }
        }

        public bool IsAlternativeStrengthUsed
        {
            get
            {
                return FIsAlternativeStrengthUsed;
            }
            set
            {
                FIsAlternativeStrengthUsed = value;
            }
        }

        public double GammaWater
        {
            get
            {
                return FGammaWater;
            }
            set
            {
                FGammaWater = value;
            }
        }

        public TWaterLineRecord[] WaterLines
        {
            get
            {
                return FWaterLines;
            }
            set
            {
                FWaterLines = value;
            }
        }

        public TPoints[][] PLLines
        {
            get
            {
                return FPLLines;
            }
            set
            {
                FPLLines = value;
            }
        }

        public TEarthquakeData EarthquakeData
        {
            get
            {
                return FEarthquakeData;
            }
            set
            {
                FEarthquakeData = value;
            }
        }

        public TUniformLoad[] UniformLoad
        {
            get
            {
                return FUniformLoad;
            }
            set
            {
                FUniformLoad = value;
            }
        }

        public TLineLoad[] LineLoad
        {
            get
            {
                return FLineLoad;
            }
            set
            {
                FLineLoad = value;
            }
        }

        public TStressTable[] StressTable
        {
            get
            {
                return FStressTable;
            }
            set
            {
                FStressTable = value;
            }
        }

        public TStabilitySideSet StabilitySide
        {
            get
            {
                return FStabilitySide;
            }
            set
            {
                FStabilitySide = value;
            }
        }

        // Degree of consolidation
        public double[] DegreeOfConsolidationLoads
        {
            get
            {
                return FDegreeOfConsolidationLoads;
            }
            set
            {
                FDegreeOfConsolidationLoads = value;
            }
        }

        public double[] DegreeOfConsolidationEarthQuake
        {
            get
            {
                return FDegreeOfConsolidationEarthQuake;
            }
            set
            {
                FDegreeOfConsolidationEarthQuake = value;
            }
        }

        public double[][] DegreeOfConsolidationSoil
        {
            get
            {
                return FDegreeOfConsolidationSoil;
            }
            set
            {
                FDegreeOfConsolidationSoil = value;
            }
        }

        public bool IsGDInternalVersion
        {
            get
            {
                return FIsGDInternalVersion;
            }
            set
            {
                FIsGDInternalVersion = value;
            }
        }

        public bool IsProbabilisticCalculation
        {
            get
            {
                return FIsProbabilisticCalculation;
            }
            set
            {
                FIsProbabilisticCalculation = value;
            }
        }

        public int LeftNumber
        {
            get
            {
                return FLeftNumber;
            }
            set
            {
                FLeftNumber = value;
            }
        }

        public int RightNumber
        {
            get
            {
                return FRightNumber;
            }
            set
            {
                FRightNumber = value;
            }
        }

        public TInterfaceCalls InterfaceCalls
        {
            get
            {
                return FInterfaceCalls;
            }
            set
            {
                FInterfaceCalls = value;
            }
        }

        public void FillSlices(int NSlices, ref Tslice[] Slices)
        {
            int i = 0;
            int StressTablePoint = 0;
            int StressTablePoint2 = 0;
            double LDeltaValue = 0;
            var LPseudoFac = new double[0];
            Tslice LSlice;
            TStabSoilData LSoil;
            for (i = 0; i < NSlices; i++)
            {
                Slices[i].addCalculatedProperties();
            }

            FNSoils = SoilData.Length;
            // Initialise pseudo factors
            LPseudoFac = new double[FNSoils];
            for (i = 0; i < FNSoils; i++)
            {
                LPseudoFac[i] = 0.0;
            }
            for (i = 0; i < NSlices; i++)
            {
                if (FIsUpliftCalculation && (i < FRightNumber) && (i > FLeftNumber))
                {
                    // in case of Uplift calculation the most unfavourable plane is searched for,
                    // a plane of 2 cm  above the slide plane is compared with a plane 2 cm below the slide plane
                    LDeltaValue = 0.02;
                    FillSliceWithData(ref Slices[i], ref StressTablePoint, LDeltaValue);
                    LSlice = new Tslice(Slices[i]);
                    LDeltaValue = -0.02;
                    FillSliceWithData(ref LSlice, ref StressTablePoint2, LDeltaValue);
                    if ((DetermineShear(Slices[i]) > DetermineShear(LSlice)))
                    {
                        Slices[i].CopySliceData(LSlice);
                        StressTablePoint = StressTablePoint2;
                    }
                }
                else
                {
                    LDeltaValue = 0.0;
                    FillSliceWithData(ref Slices[i], ref StressTablePoint, LDeltaValue);
                }
                if (Slices[i].StrengthType == TMatStrengthTypeSet.mstStressTab || Slices[i].StrengthType == TMatStrengthTypeSet.mstPseudoStressTab)
                {
                    LSoil = SoilData[Slices[i].SoilNumber];
                    // voor de pseudofactor wordt de lengte in het materiaal bijgehouden
                    LPseudoFac[Slices[i].SoilNumber] = LPseudoFac[Slices[i].SoilNumber] + Slices[i].ArcLength*
                                                       LSoil.StressTable.Delta(Slices[i].EffectiveStress, StressTablePoint);
                }
            }
            // nu per lamel voor de pseudo de coh en phi bepalen
            if ((FDefaultStrength == TMatStrengthTypeSet.mstPseudoStressTab))
            {
                Det_Tau_pseudo(ref LPseudoFac, ref Slices, NSlices);
            }
        }

        // (*
        // {=======================================================================================================================
        //
        // Date        ID   Modification
        // 2004-02-19  Best Created
        // =======================================================================================================================}
        // procedure TSliceFill.InitialiseStresstablefactor;
        // var
        // i: Integer;
        // LStandard: Boolean;
        // LInputType: TProbInputValuesSet;
        // LCalculationType: TProbCalculationType;
        // LUseProbDefaults: Boolean;
        // LConvert: TConvertTypeSet;
        // LSoil: TStabSoilData;
        // LSoilStrengthType: TMatStrengthTypeSet;
        // begin
        // for i := Low(FSoilData) to High(FSoilData) do
        // begin
        // LSoil := FSoilData[i];
        // LSoil.StressTableFactor := CNormValue;
        // end;
        // end;
        //
        // if (FIsProbabilisticCalculation) then
        // begin
        // LInputType := LSoil.ProbInputValues;
        // if (LInputType = pivNone) then
        // begin
        // LInputType := FDefaultInputValues;
        // end;
        // LUseProbDefaults := LSoil.UseProbDefaults;
        // LCalculationType := FProbabilisticCalculationType;
        //
        // LSoilStrengthType := LSoil.MatStrengthType;
        // if (LSoilStrengthType = mstNone) then
        // begin
        // LSoilStrengthType := GetStrengthType;
        // end;
        //
        // if (LSoilStrengthType = mstStressTab) then
        // begin
        //
        // {... calculation when only one columns is given ...}
        // if (LUseProbDefaults) then
        // begin
        // LStandard := (LCalculationType = pctDeterministicMean) and (LInputType =
        // pivMean);
        // if (not LStandard) then
        // begin
        // LStandard := (LCalculationType = pctDeterministicDesign) and
        // (LInputType = pivDesign);
        // end;
        // if (not LStandard) then
        // begin
        // LStandard := (LCalculationType = pctProbabilistic) and (LInputType =
        // pivMean);
        // end;
        //
        // {Get the factor for calculation one into the other}
        // if (not LStandard) then
        // begin
        // LConvert := ctStressTabFactor;
        // LSoil.StressTableFactor := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
        // end;
        // end;
        //
        // end;
        // end;
        // end;
        //
        // end;
        // *)
        // =======================================================================================================================
        // Date        ID   Modification
        // 2008-04-23  Best Created
        // =======================================================================================================================
        private TMatStrengthTypeSet GetStrengthType()
        {
            TMatStrengthTypeSet result;
            result = TMatStrengthTypeSet.mstCphi;
            switch (FDefaultStrength)
            {
                case TMatStrengthTypeSet.mstCphi:
                    result = TMatStrengthTypeSet.mstCphi;
                    break;
                case TMatStrengthTypeSet.mstStressTab:
                    result = TMatStrengthTypeSet.mstStressTab;
                    break;
                case TMatStrengthTypeSet.mstPseudoStressTab:
                    result = TMatStrengthTypeSet.mstPseudoStressTab;
                    break;
                case TMatStrengthTypeSet.mstMeasuredCu:
                    result = TMatStrengthTypeSet.mstMeasuredCu;
                    break;
                case TMatStrengthTypeSet.mstGradientCu:
                    result = TMatStrengthTypeSet.mstGradientCu;
                    break;
                case TMatStrengthTypeSet.mstCalculatedCu:
                    result = TMatStrengthTypeSet.mstCalculatedCu;
                    break;
                case TMatStrengthTypeSet.mstNone:
                    result = TMatStrengthTypeSet.mstCphi;
                    break;
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2008-06-05  Best Created
        // =======================================================================================================================
        private double GetCuValues(ref Tslice ASlice)
        {
            double result;
            double LCuGem, LCuActief, LCuPasief;
            // LCuActiveTop, LCuPassiveTop: Double;
            // LCuActiveBottom, LCuPassiveBottom: Double;
            // LCuTop, LCuBottom: Double;
            double LPop;
            TStabSoilData LSoil;
            // LConvert: TConvertTypeSet;
            double LUsedStress;
            double LZLayerTop = 0;
            double LZLayerBottom = 0;
            result = 0;
            LSoil = SoilData[ASlice.SoilNumber];
            ASlice.IsInActiveArea = ((FStabilitySide == TStabilitySideSet.RightSide) && (ASlice.ZBottomRight < ASlice.ZBottomLeft)) || ((FStabilitySide == TStabilitySideSet.LeftSide) && (ASlice.ZBottomRight > ASlice.ZBottomLeft));
            LPop = LSoil.POP;
            // ... Cu Measured ...
            if ((ASlice.StrengthType == TMatStrengthTypeSet.mstMeasuredCu))
            {
                ASlice.Vertical.DetermineTopAndBottomFromLayer(ASlice.ZBottomMid, ref LZLayerTop, ref LZLayerBottom);
                if (LZLayerTop > LZLayerBottom)
                {
                    ASlice.InterpolationFactorCu = (ASlice.ZBottomMid - LZLayerBottom)/(LZLayerTop - LZLayerBottom);
                }
                else
                {
                    ASlice.InterpolationFactorCu = 1;
                }
                result = ASlice.InterpolationFactorCu*(LSoil.CuTop - LSoil.CuBottom) + LSoil.CuBottom;

                if (FIsProbabilisticCalculation)
                {
                    LCuGem = ASlice.InterpolationFactorCu*(LSoil.CuTop - LSoil.CuBottom) + LSoil.CuBottom;
                    LCuActief = ASlice.InterpolationFactorCu*(LSoil.CuActiveTop - LSoil.CuActiveBottom) + LSoil.CuActiveBottom;
                    LCuPasief = ASlice.InterpolationFactorCu*(LSoil.CuPassiveTop - LSoil.CuPassiveBottom) + LSoil.CuPassiveBottom;

                    if (ASlice.IsInActiveArea)
                    {
                        result = (LCuGem + (LCuActief - LCuGem)*Math.Sin(2*Math.Abs(ASlice.BottomAngle)));
                    }
                    else
                    {
                        // passive area
                        result = (LCuGem - (LCuGem - LCuPasief)*Math.Sin(2*Math.Abs(ASlice.BottomAngle)));
                    }
                }
                else
                {
                    // {... not probabilistic ...}
                    result = ASlice.InterpolationFactorCu*(LSoil.CuTop - LSoil.CuBottom) + LSoil.CuBottom;
                }
            }
            // begin
            // LConvert := ctCuActiveTop;
            // LCuActiveTop := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
            // LConvert := ctCuActiveBottom;
            // LCuActiveBottom := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
            //
            // LConvert := ctCuPassiveTop;
            // LCuPassiveTop := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
            // LConvert := ctCuPassiveBottom;
            // LCuPassiveBottom := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
            //
            // LConvert := ctCuTop;
            // LCuTop := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
            // LConvert := ctCuBottom;
            // LCuBottom := GetSoilParameterValue(LConvert, LSoil, MyGlobals);
            //
            // LCuGem := InterpolationFactorCu * (LCuTop - LCuBottom) + LCuBottom;
            // LCuActief := InterpolationFactorCu * (LCuActiveTop - LCuActiveBottom) +
            // LCuActiveBottom;
            // LCUPasief := InterpolationFactorCu * (LCuPassiveTop - LCuPassiveBottom) +
            // LCuPassiveBottom;
            //
            // if (IsInActiveArea) then
            // Result := (LCuGem + (LCuActief - LCuGem) * Sin(2 * Abs(BottomAngle)))
            // else
            // Result := (LCuGem - (LCuGem - LCuPasief) * Sin(2 * Abs(BottomAngle)))
            // end
            // else
            // {... not probabilistic ...}
            // begin
            // Result := InterpolationFactorCu * (CuTop - CuBottom) + CuBottom;
            // end;
            // end;
            // *)

            // ... Cu Gradient (not yet probabilistic) ...
            if ((ASlice.StrengthType == TMatStrengthTypeSet.mstGradientCu))
            {
                ASlice.Vertical.DetermineTopAndBottomFromLayer(ASlice.ZBottomMid, ref LZLayerTop, ref LZLayerBottom);
                result = LSoil.CuTop + (LZLayerTop - ASlice.ZBottomMid)*LSoil.CuGradient;
            }
            // ... Cu calculated ...
            if ((ASlice.StrengthType == TMatStrengthTypeSet.mstCalculatedCu))
            {
                LUsedStress = Math.Max(ASlice.PreLoadStress + LPop, ASlice.EffectiveStress);
                if (FIsProbabilisticCalculation)
                {
                    if (ASlice.IsInActiveArea)
                    {
                        result = (LUsedStress)*(LSoil.RatioCuPc + (LSoil.RatioCuPcActive - LSoil.RatioCuPc)*
                                                Math.Sin(2*Math.Abs(ASlice.BottomAngle)));
                    }

                    else
                    {
                        result = (LUsedStress)*(LSoil.RatioCuPc - (LSoil.RatioCuPc - LSoil.RatioCuPcPassive)*
                                                Math.Sin(2*Math.Abs(ASlice.BottomAngle)));
                    }
                }
                else
                {
                    // ... not probabilistic ...
                    result = LUsedStress*LSoil.RatioCuPc;
                }
            }
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2008-04-23  Best Created
        // =======================================================================================================================
        private void FillSliceWithData(ref Tslice ASlice, ref int StressTablePoint, double DeltaValue)
        {
            const double CAlmostZero = 0.001;
            TStabSoilData LSoil;
            double LForce = 0;
            double LStressAtBot = 0;
            double LTotUniformLoadStress = 0;
            // Stress due to uniform load
            double LUnifLoadstress = 0;
            double LPointStress = 0;
            // spanning tgv puntlast(en)
            double Laanpas = 0;
            double LTopZ = 0;
            double LCohesion = 0;
            double LPhi = 0;
            TVertical LVertical;
            double LDOCQuake = 0;
            double LZBottom = 0;
            ASlice.Vertical = new TVertical();
            LVertical = ASlice.Vertical;
            LZBottom = ASlice.ZBottomMid + DeltaValue;
            LVertical.SoilData.AddRange(FSoilData);
            LVertical.XCoordinate = ASlice.XMid;
            LVertical.ZCoordinate = ASlice.ZBottomMid;
            LVertical.PLLines = FPLLines;
            LVertical.GammaWater = FGammaWater;
            LVertical.WaterLines.AddRange(FWaterLines);
            LVertical.FreaticLine = FFreaticLine;
            LVertical.SurfaceLine = FSurfaceLine;
            LVertical.PreconStressLine = FPreconStressLine;
            LVertical.FillVerticalWithWaterData();
            LVertical.DegreeOfConsolidationSoil = FDegreeOfConsolidationSoil;
            LVertical.InterfaceCalls = FInterfaceCalls;
            ASlice.InterfaceCalls = FInterfaceCalls;
            LVertical.FillVertical();
            if ((FreaticLine != null))
            {
                LVertical.Freaticlevel = MStabDatafunctions.ZTopAtSurface(ASlice.XMid, FreaticLine);
            }
            else
            {
                LVertical.Freaticlevel = -999;
            }
            LVertical.SoilLevel = MStabDatafunctions.ZTopAtSurface(ASlice.XMid, FSurfaceLine);
            ASlice.SoilNumber = ASlice.InterfaceCalls.DetermineSoilNumberFromPoint(ASlice.XMid, LZBottom);
            ASlice.InLayer = ASlice.InterfaceCalls.DetermineLayerNumberFromPoint(ASlice.XMid, LZBottom);
            ASlice.LayersAbove = LVertical.DetermineSoilsAbove(LZBottom, ASlice.InLayer);
            LSoil = SoilData[ASlice.SoilNumber];

            ASlice.SoilGroupNumber = LSoil.SoilGroup;
            // pseudostress only possible if all mateials have pseudo stress tables
            // so material is overruled
            if ((FDefaultStrength == TMatStrengthTypeSet.mstPseudoStressTab))
            {
                ASlice.StrengthType = TMatStrengthTypeSet.mstPseudoStressTab;
            }
            else
            {
                // get the Material strength type (default if StrengthType = mstNone)
                ASlice.StrengthType = LSoil.MaterialStrength;
                if (
                    (new ArrayList(new TMatStrengthTypeSet[]
                        {
                            TMatStrengthTypeSet.mstNone, TMatStrengthTypeSet.Default
                        })
                        .Contains(ASlice.StrengthType)))
                {
                    ASlice.StrengthType = GetStrengthType();
                }
            }
            if (
                (new ArrayList(new TMatStrengthTypeSet[]
                    {
                        TMatStrengthTypeSet.mstStressTab, TMatStrengthTypeSet.mstPseudoStressTab
                    }).Contains(
                        ASlice.StrengthType)))
            {
                ASlice.SpanTabel = LSoil.StressTableNr;
            }
            else
            {
                ASlice.SpanTabel = 0;
            }

            ASlice.ZFrea = LVertical.Freaticlevel;

            // Calculate pore on surface
            ASlice.PoreOnSurface = LVertical.DeterminePoreOnSurface();
            try
            {
                // Calculate total water force on slice
                LForce = ASlice.PoreOnSurface*ASlice.Width/Math.Cos(ASlice.TopAngle);
            }
            catch
            {
                LForce = 0;
            }
            // HPoreOnSurface  = horizontale waterkracht op maaiveld
            ASlice.HPoreOnSurface = LForce*Math.Sin(ASlice.TopAngle);
            // VPoreOnSurface  = vertikale waterkracht op maaiveld
            ASlice.VPoreOnSurface = -LForce*Math.Cos(ASlice.TopAngle);
            // gewicht lamel = grondspan * breedte
            // LStressAtBot := FInterfacecalls.DetermineTotalstress(XMid, ZBottomMid);
            LStressAtBot = LVertical.DetermineTotalstress(ASlice.ZBottomMid);
            ASlice.Weight = (LStressAtBot - ASlice.PoreOnSurface)*ASlice.Width;
            // ... Influence of earth quake ...
            // Find stresse due to quake
            ASlice.SigmaSoilQuake = 0;
            ASlice.SigmaWaterQuake = 0;
            if ((Math.Abs(FEarthquakeData.VerticalQuakeFactor) > 0.00001))
            {
                ASlice.SigmaSoilQuake = (LStressAtBot - ASlice.PoreOnSurface)*FEarthquakeData.VerticalQuakeFactor;
                LDOCQuake = EarthquakeData.DegreeOfConsolidation[ASlice.InLayer];
                ASlice.SigmaWaterQuake = LVertical.DeterminePoreFromVertQuake(ASlice.ZBottomMid, ASlice.SigmaSoilQuake,
                                                                              LDOCQuake,
                                                                              EarthquakeData.VerticalQuakeFactor);
            }
            // find horizontal force due to quake
            ASlice.HorizontalQuakeforce = 0;
            if ((FEarthquakeData.HorizontalQuakeFactor > 0.00001))
            {
                ASlice.HorizontalQuakeforce = ASlice.Weight*FEarthquakeData.HorizontalQuakeFactor;
            }
            // ... Stresses due to uniform loads ...
            LTotUniformLoadStress = 0;
            ASlice.WaterLoadStress = 0;
            for (int i = 0; i < FUniformLoad.Length; i++)
            {
                TUniformLoad _wvar1 = FUniformLoad[i];
                LUnifLoadstress = _wvar1.GetStressAtPosition(ASlice.XMid, ASlice.ZBottomMid);
                // each slice can have one or more external loads
                // find excess pore pressure due to tmp loads
                if ((LUnifLoadstress > 0))
                {
                    if ((_wvar1.LoadType == TLoadTypeSet.Temporary))
                    {
                        if ((ASlice.ZBottomMid < ASlice.ZFrea) || FIsCapillairLoadUsed)
                        {
                            Laanpas = DegreeOfConsolidationLoads[ASlice.InLayer]/100;
                            ASlice.WaterLoadStress = ASlice.WaterLoadStress + (1 - Laanpas)*LUnifLoadstress;
                        }
                    }
                }
                LTotUniformLoadStress = LTotUniformLoadStress + LUnifLoadstress;
            }
            ASlice.ExternalLoad = LTotUniformLoadStress*ASlice.Width;
            // ... Stresses due to point-loads ...
            LPointStress = 0.0;
            for (int i = 0; i < ASlice.LineLoad.Count; i++)
            {
                TLineLoad _wvar2 = ASlice.LineLoad[i];
                LPointStress = LPointStress + _wvar2.GetStressAtLamel(ASlice.XMid, ASlice.ZBottomMid, ASlice.Width);
            }
            // ... Add all stresses ...
            ASlice.TotalStress = LStressAtBot + LTotUniformLoadStress + LPointStress + ASlice.SigmaSoilQuake;
            // due to negative uniform load stress total stress can become negative
            if ((ASlice.TotalStress < 0))
            {
                // add an extra negative pore pressure component ?
                ASlice.TotalStress = 0;
            }

            ASlice.AanpasPore = LVertical.DetermineDegreeOfConsolidationSoilStress(ASlice.ZBottomMid, ASlice.InLayer);
            // FInterfaceCalls.DeterminePorePressureDegreeOfConsolidationLayers(XMid, ZBottomMid);
            // PnLinePore := FInterfaceCalls.DeterminePorePressurePiezometricLevels(XMid, ZBottomMid);
            ASlice.PnLinePore = LVertical.DetermineWaternetStress(ASlice.ZBottomMid);
            ASlice.TotalPore = ASlice.PnLinePore + ASlice.AanpasPore + ASlice.SigmaWaterQuake + ASlice.WaterLoadStress;
            // Internal a research may be performed
            if (FIsGDInternalVersion)
            {
                ASlice.TotalPore = LSoil.PorePressureFactor*(ASlice.TotalPore + LSoil.ExcessPorePressure);
            }

            ASlice.EffectiveStress = ASlice.TotalStress - ASlice.TotalPore;
            if (ASlice.EffectiveStress < 0)
            {
                ASlice.TotalPore = ASlice.TotalPore + ASlice.EffectiveStress;
                ASlice.EffectiveStress = 0.0;
            }
            if (ASlice.ZFrea > ASlice.ZBottomMid)
            {
                ASlice.HydrostaticPore = (ASlice.ZFrea - ASlice.ZBottomMid)*FGammaWater;
            }
            else
            {
                ASlice.HydrostaticPore = 0;
            }

            // alleen nodig voor aardbevingen
            ASlice.Zwaartepunt = LVertical.DeterminCentreOfGravity(ASlice.ZBottomMid);

            if (
                new ArrayList(new TMatStrengthTypeSet[]
                    {
                        TMatStrengthTypeSet.mstCalculatedCu
                    }).Contains(
                        ASlice.StrengthType))
            {
                // ... Get (initial)stresses for cu calculated ...
                LTopZ = ASlice.InterfaceCalls.DetermineZReferencelevelCuCalculated(ASlice.XMid);
                if (Math.Abs(LTopZ - ASlice.ZTopMid) > CAlmostZero)
                {
                    if (LTopZ > ASlice.ZBottomMid)
                    {
                        ASlice.PreLoadStress = LVertical.DeterminePreLoadStress(ASlice.ZBottomMid, LTopZ) -
                                               ASlice.PnLinePore;
                        // PreLoadStress := Sigma_totIncrement(XMid, LTopY, YBottomMid) - aSlice.PnLinePore;
                        // ASlice.PreLoadStress = FInterfaceCalls.DeterminePreConsolidationTotalStress(XMid, ZBottomMid, LTopZ) -
                        // ASlice.PnLinePore;
                        if ((ASlice.PreLoadStress < 0))
                        {
                            ASlice.PreLoadStress = 0;
                        }
                    }
                    else
                    {
                        ASlice.PreLoadStress = LStressAtBot - ASlice.PnLinePore;
                    }
                }
                else
                {
                    ASlice.PreLoadStress = LStressAtBot - ASlice.PnLinePore;
                }
            }

            // de coh en phi
            // Cu waarde aan onderkant lamel bij Cu meaured
            // only valid if EffectiveStress is calculated
            if (
                new ArrayList(new TMatStrengthTypeSet[]
                    {
                        TMatStrengthTypeSet.mstMeasuredCu, TMatStrengthTypeSet.mstCalculatedCu,
                        TMatStrengthTypeSet.mstGradientCu
                    }).Contains(ASlice.StrengthType))
            {
                ASlice.Cu = GetCuValues(ref ASlice);
                ASlice.PhiBottom = 0;
                ASlice.CohBottom = ASlice.Cu;
            }
            else if (
                (!(new ArrayList(new TMatStrengthTypeSet[]
                    {
                        TMatStrengthTypeSet.mstStressTab, TMatStrengthTypeSet.mstPseudoStressTab
                    })
                      .Contains(ASlice.StrengthType))))
            {
                ASlice.PhiBottom = LSoil.Phi;
                ASlice.CohBottom = LSoil.Cohesion;
            }
            // Initialiseer Sigma-alfa met de vertikale grondspanning
            if (
                new ArrayList(new TMatStrengthTypeSet[]
                    {
                        TMatStrengthTypeSet.mstStressTab, TMatStrengthTypeSet.mstPseudoStressTab
                    }).
                    Contains(ASlice.StrengthType))
            {
                ASlice.SigmaAlfa = ASlice.EffectiveStress;
                // cohesie en fie aan onderkant lamel
                LSoil.GetCohAndPhiFromTable(ASlice.SigmaAlfa, ASlice.PseudoFactor, ref LCohesion, ref LPhi,
                                            ref StressTablePoint);
                ASlice.CohBottom = LCohesion;
                ASlice.PhiBottom = LPhi;
            }
            if ((ASlice.StrengthType == TMatStrengthTypeSet.mstCphi) && FIsAlternativeStrengthUsed)
            {
                ASlice.Dilatancy = LSoil.Dilatancy;
            }
            else
            {
                ASlice.Dilatancy = ASlice.PhiBottom;
            }
        }

        // =======================================================================================================================
        // Description: The ksi from each layer is calculated
        // For background see
        // Pseudo karakteristiek waarden
        // (aanpassing t.b.v. Windows versie van MStab)
        // DF 510090.01  (Calle)
        //
        // Date        ID   Modification
        // 2002-02-21  Best Created
        // =======================================================================================================================
        private void CalculatePseudofactors(ref double[] APseudoFac)
        {
            int RowIndex;
            int i;
            int j;
            int LSoilType;
            TStressTable LSpanTabel;
            TStressTable LInitSpanTabel;
            int LNMeasure;
            int LOtherMeasure;
            double LCorrelation;
            int LNRows;
            int LNColumns;
            TMatrix LCorrelationMatrix;
            TMatrix LStdvVector;
            TMatrix LStdvVectorTransposed;
            TMatrix LTempVector;
            TMatrix LNoemerMatrix;
            double LNoemerResult;
            // ... initialise matrices ...
            LStdvVector = null;
            LTempVector = null;
            LNoemerMatrix = null;
            LStdvVectorTransposed = null;
            // ... define  correlation matrix ...
            LNRows = FNSoils;
            LNColumns = FNSoils;
            LCorrelationMatrix = new TMatrix();
            LCorrelationMatrix.SetDimensions(LNRows, LNColumns);
            LCorrelationMatrix.SetToNaught();
            // ... fill correlation Matrix ...
            for (i = 1; i <= LNColumns; i++)
            {
                LSoilType = SoilData[i - 1].SoilGroup;
                LInitSpanTabel = SoilData[i - 1].StressTable;
                LNMeasure = LInitSpanTabel.NTests;
                for (j = i; j <= LNRows; j++)
                {
                    if ((i == j))
                    {
                        LCorrelationMatrix[j - 1, i - 1] = 1;
                    }
                    else
                    {
                        if ((LSoilType > 0))
                        {
                            if ((SoilData[j - 1].SoilGroup == LSoilType))
                            {
                                LSpanTabel = SoilData[j - 1].StressTable;
                                // ... if stresstables are the same correlation is 1 ...
                                if ((LInitSpanTabel == LSpanTabel))
                                {
                                    LCorrelation = 1;
                                }
                                else
                                {
                                    LOtherMeasure = LSpanTabel.NTests;
                                    LCorrelation = 1/Math.Sqrt((1 + 1/(0.25*LNMeasure))*(1 + 1/(0.25*LOtherMeasure)));
                                }
                                // ... symetrische matrix ...
                                LCorrelationMatrix[i - 1, j - 1] = LCorrelation;
                                LCorrelationMatrix[j - 1, i - 1] = LCorrelation;
                            }
                        }
                    }
                }
            }
            // ... create the standard deviation vector ...
            LNRows = FNSoils;
            LNColumns = 1;
            LStdvVector = new TMatrix();
            LStdvVector.SetDimensions(LNRows, LNColumns);
            LStdvVector.SetToNaught();
            // ... fill the standard deviation vector ...
            for (RowIndex = 0; RowIndex < LNRows; RowIndex++)
            {
                LStdvVector[RowIndex, LNColumns - 1] = APseudoFac[RowIndex];
            }
            // ... create the tranposed standard deviation vector ...
            LNRows = 1;
            LNColumns = FNSoils;
            LStdvVectorTransposed = new TMatrix();
            LStdvVectorTransposed.SetDimensions(LNRows, LNColumns);
            LStdvVectorTransposed.SetToNaught();
            LStdvVectorTransposed = LStdvVector.Transpose();
            // ... create a temporary vector ...
            LTempVector = LCorrelationMatrix.Multiply(LStdvVector);
            LNoemerMatrix = LStdvVectorTransposed.Multiply(LTempVector);
            LNoemerResult = Math.Sqrt(Convert.ToDouble(LNoemerMatrix[0, 0]));
            LNColumns = LTempVector.ColumnCount;
            Debug.Assert(LNColumns == 1);
            for (i = 0; i < FNSoils; i++)
            {
                if ((LNoemerResult == 0))
                {
                    APseudoFac[i] = 0;
                }
                else
                {
                    APseudoFac[i] = LTempVector[i, LNColumns - 1]/LNoemerResult;
                }
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2001-05-28  Best Created
        // =======================================================================================================================
        private void GetLocalPseudoFactors(ref double[] APseudoFac)
        {
            int i;
            double T;
            // ...  Normaliseer de PseudoFactor's ...
            T = 0.0;
            for (i = 0; i < FNSoils; i++)
            {
                T = T + Math.Pow(APseudoFac[i], 2);
            }
            if (T <= 0)
            {
                T = 0.0;
            }
            else
            {
                T = Math.Sqrt(T);
            }
            if (T > 1.0E-30)
            {
                for (i = 0; i < FNSoils; i++)
                {
                    APseudoFac[i] = APseudoFac[i]/T;
                }
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2008-06-11  Best Created
        // =======================================================================================================================
        private void Det_Tau_pseudo(ref double[] APseudoFac, ref Tslice[] Slices, int NSlices)
        {
            int i = 0;
            int TsPoint = 0;
            double kor = 0;
            double coh = 0;
            double fie = 0;
            double t = 0;
            TStabSoilData LSoil;

            // hulp parameters
            if (FIsGlobalMesurement)
            {
                CalculatePseudofactors(ref APseudoFac);
            }
            else
            {
                GetLocalPseudoFactors(ref APseudoFac);
            }
            // Voor elke lamel, bepaal de lamelgrootheden
            for (i = 0; i < NSlices; i++)
            {
                LSoil = SoilData[Slices[i].SoilNumber];

                // lamel^[i,15]  = korrelspanning aan onderkant lamel
                kor = Slices[i].EffectiveStress;
                // #best stond er        PseudoFactor := PseudoFac[InLayer];
                Slices[i].PseudoFactor = APseudoFac[Slices[i].SoilNumber];
                // Bepaal Pseudo karakteristiek
                LSoil.StressTable.GetCohPhi(kor, Slices[i].PseudoFactor, ref coh, ref fie, ref TsPoint);
                // cohesie aan onderkant lamel
                Slices[i].CohBottom = coh;
                // fie aan onderkant lamel (in graden)
                Slices[i].PhiBottom = fie;
                // Tau langs onderkant lamel
                t = fie*Math.PI/180.0;
                t = Math.Sin(t)/Math.Cos(t);
                Slices[i].NormalStress = kor;
                Slices[i].ShearStress = kor*t + coh;
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2008-08-11  Best Created
        // =======================================================================================================================
        private double DetermineShear(Tslice aSlice)
        {
            double result;
            result = aSlice.EffectiveStress*aSlice.PhiBottom*Math.PI/180.0 + aSlice.CohBottom;
            return result;
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2008-04-23  Best Created
        // =======================================================================================================================
    }

    // end TSliceFill
}