using System;
using System.Collections.Generic;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Consolidation;
using Deltares.Geotechnics.Soils;
using Deltares.Standard.Language;

//
//  Interpolation columns consists of x left and x right
//  There is column data and area data
//  column data is:   Level freatic line
//                    Load on top (water load and other load)
//                    Top angle of the slice  
//                  
// 
//   area data is data needed at the bottom of the slice
// 
//  
//  loads on surface  this includes the uniform loads 
//  Depending on the geometry a column has 0ne or more areas
//  Each area has 4 corners in which the data (stresses) is known
//  
// 
//  calculation of the data
//  as the columns are between all the points in the geometry the data is 
//  best calculated next to the column boundarys
// 
//  This is done at 1/3 and 2/3 of the column width
//  by extrapolationdata is found in the corners of the areas 
//  cq boundaries of the columns
// 
//

namespace Deltares.Stability.Calculation2
{
    public class InterpolationColumn
    {
        private const double delta = 0.00001;
        private const double CDegreeToRad = Math.PI/180.0;
        private const double CRadToDegree = 180.0/Math.PI;

        //  ... column data

        private List<InterpolationArea> areas = new List<InterpolationArea>();

        // help variables
        private double helpLeftX;
        private List<double> helpLeftpoints;
        private double helpPoreLeft = 0.0;
        private double helpPoreRight = 0.0;
        private SoilProfile1D helpProfileLeft;
        private SoilProfile1D helpProfileRight;
        private double helpRightX;
        private List<double> helpRightpoints;

        private double helpZFreaLeft;
        private double helpZFreaRight;
        private double helpZPreconBounLeft = double.NaN;
        private double helpZPreconBounRight = double.NaN;
        private double xCoordinate;
        private double zCoordinate;
        public InterpolationColumn() {}

        public InterpolationColumn(InterpolationColumn aColumn)
        {
            for (int i = 0; i < aColumn.areas.Count; i++)
            {
                var area = new InterpolationArea(aColumn.areas[i]);
                areas.Add(area);
            }
        }

        public double TopLoadPermLeft { get; set; }
        public double TopLoadPermRight { get; set; }
        public double TopLoadTempLeft { get; set; }
        public double TopLoadTempRight { get; set; }
        public double PoreOnSurfaceLeft { get; set; }
        public double PoreOnSurfaceRight { get; set; }
        public double ZFreaLeft { get; set; }
        public double ZFreaRight { get; set; }
        public double ZPreConBoundaryLeft { get; set; }
        public double ZPreConBoundaryRight { get; set; }

        public double XLeft { get; set; }
        public double XRight { get; set; }
        public double GammaWater { get; set; }

        //  ... geometry data from slice and horizontal interpolation factor ..
        //

        public double HorInpolFactor { get; set; }
        //
        //  ... calculated column data ...
        //  these values are calculated by (horizontal) interpolation 
        //  in the method XIsInColumn
        //
        public double UniformLoadPermMagnitude { get; set; }
        public double UniformLoadTempMagnitude { get; set; }
        public double ZFrea { get; set; }
        public double TopAngle { get; set; }
        public double PoreOnSurface { get; set; }
        //     public double HPoreOnSurface { get; set; }
        //    public double VPoreOnSurface { get; set; }

        public List<InterpolationArea> Areas
        {
            get
            {
                return areas;
            }
            set
            {
                areas = value;
            }
        }

        public bool XIsInColumn(double AXCoor)
        {
            if ((AXCoor >= XLeft) && (AXCoor <= XRight))
            {
                //
                // column found calculate the horizontal interpolation factor
                // and some column parameters as well
                //
                xCoordinate = AXCoor;
                HorInpolFactor = (xCoordinate - XLeft)/(XRight - XLeft);
                return true;
            }
            return false;
        }

        public void DetermineColumnData()
        {
            if (Areas.Count > 0)
            {
                double zLeft = Areas[0].TopLeftPoint.ZCoordinate;
                double zRight = Areas[0].TopRightPoint.ZCoordinate;
                TopAngle = Math.Atan2((zRight - zLeft), (XRight - XLeft));
            }
        }

        public void RetrieveColumnData()
        {
            PoreOnSurface = PoreOnSurfaceLeft + HorInpolFactor * (PoreOnSurfaceRight - PoreOnSurfaceLeft);
            //        HPoreOnSurface = PoreOnSurface * Math.Sin(TopAngle) ;
            //        VPoreOnSurface = -PoreOnSurface * Math.Cos(TopAngle);
            ZFrea = ZFreaLeft + HorInpolFactor*(ZFreaRight - ZFreaLeft);
            UniformLoadPermMagnitude = TopLoadPermLeft + HorInpolFactor * (TopLoadPermRight - TopLoadPermLeft);
            UniformLoadTempMagnitude = TopLoadTempLeft + HorInpolFactor * (TopLoadTempRight - TopLoadTempLeft);
        }

        /// <summary>
        ///  Get the area with th x and y coordinate
        /// </summary>
        public InterpolationArea getArea(double aZCoordinate)
        {
            zCoordinate = aZCoordinate;
            double ZTop;
            double ZBottom;
            // Find the array for the ZValue 
            foreach (var LocalArea in Areas)
            {
                ZTop = LocalArea.TopLeftPoint.ZCoordinate + HorInpolFactor*
                       (LocalArea.TopRightPoint.ZCoordinate - LocalArea.TopLeftPoint.ZCoordinate);
                ZBottom = LocalArea.BottomLeftPoint.ZCoordinate + HorInpolFactor*
                          (LocalArea.BottomRightPoint.ZCoordinate - LocalArea.BottomLeftPoint.ZCoordinate);

                if ((zCoordinate <= ZTop + delta) && (zCoordinate >= ZBottom - delta))
                {
                    //  ... area is found add the interpolation factors to this area ...
                    //
                    LocalArea.HorizontalInterpolationFactor = HorInpolFactor;
                    LocalArea.VerticalInterpolationFactor = (zCoordinate - ZTop)/(ZBottom - ZTop);
                    LocalArea.ZCoordinate = zCoordinate;
                    return LocalArea;
                }
            }
            return null;
        }

        /// <summary>
        /// Get the index (starting at 0) of the area at a certain Z coordinate
        /// </summary>
        /// <param name="aZCoordinate">Z coordinate</param>
        /// <returns>The index of the area</returns>
        public int DetermineAreaNumber(double aZCoordinate)
        {
            int i;
            double ZTop;
            double ZBottom;

            // Find the array for the ZValue 
            for (i = 0; i < Areas.Count; i++)
            {
                ZTop = Areas[i].TopLeftPoint.ZCoordinate + HorInpolFactor*(Areas[i].TopRightPoint.ZCoordinate - Areas[i].TopLeftPoint.ZCoordinate);
                ZBottom = Areas[i].BottomLeftPoint.ZCoordinate + HorInpolFactor*(Areas[i].BottomRightPoint.ZCoordinate - Areas[i].BottomLeftPoint.ZCoordinate);

                if ((aZCoordinate <= ZTop) && (aZCoordinate >= ZBottom))
                {
                    return i;
                }
            }
            return 0;
        }

        /// <summary>
        /// help verticals are created 
        /// these are private parts as they are needed everywhere in column
        /// </summary>
        /// <param name="profile2D"></param>
        public void CreateHelpverticals(SoilProfile2D profile2D)
            //========================================================
        {
            helpLeftX = RoundThreeDigits(GetXAtOneThird());
            helpRightX = RoundThreeDigits(GetXAtTwoThird());

            helpProfileLeft = profile2D.GetSoilProfile1D(helpLeftX);
            helpProfileRight = profile2D.GetSoilProfile1D(helpRightX);
        }

        public void FillPointsWithData(Waternet waternet)
            //================================================
        {
            RetrieveColumnData();
            FillPointsWithTotalStressData();
            if (waternet != null && waternet.PhreaticLine != null)
            {
                FillPointsWithWaternetPressures(waternet);
            }
            FillPointsWithPorePressuresDueToDegreeOfConsolidationLoad();
            FillPointsWithTotalPorePressures();
            FillPointsWithEffectiveStresses();
            FillPointsWithPreLoadStressData();
        }

        public void AddColumnData(Waternet waternet, GeometryPointString preloadSurfaceLine, List<UniformLoad> uniformLoads)
        {
            //TODO should waternet.PhreaticLine check in validator / default values etc?
            if (waternet != null && waternet.PhreaticLine != null)
            {
                GammaWater = waternet.UnitWeight;
                helpZFreaLeft = waternet.PhreaticLine.GetZAtX(helpLeftX);
                helpZFreaRight = waternet.PhreaticLine.GetZAtX(helpRightX);

                // DetermineAreaNumber freatic line in column
                AreaBoundary freaBoundary = GetAreaBoundaryValuesByExtrapolation(helpZFreaLeft, helpZFreaRight);
                ZFreaLeft = freaBoundary.LeftValue;
                ZFreaRight = freaBoundary.RightValue;
                if (helpZFreaLeft > helpProfileLeft.Layers[0].TopLevel)
                {
                    helpPoreLeft = (helpZFreaLeft - helpProfileLeft.Layers[0].TopLevel) * GammaWater;
                }

                if (helpZFreaRight > helpProfileRight.Layers[0].TopLevel)
                {
                    helpPoreRight = (helpZFreaRight - helpProfileRight.Layers[0].TopLevel) * GammaWater;
                }
                AreaBoundary topBoundary = GetAreaBoundaryValuesByExtrapolation(helpPoreLeft, helpPoreRight);
                PoreOnSurfaceLeft = topBoundary.LeftValue;
                PoreOnSurfaceRight = topBoundary.RightValue;
            }
            else
            {
                PoreOnSurfaceLeft = 0;
                PoreOnSurfaceRight = 0;
            }

            // determine the pre load boundaries
            if (preloadSurfaceLine != null)
            {
                helpZPreconBounLeft = preloadSurfaceLine.GetZAtX(helpLeftX);
                helpZPreconBounRight = preloadSurfaceLine.GetZAtX(helpRightX);
                // DetermineAreaNumber freatic line in column
                AreaBoundary PreConBoundary = GetAreaBoundaryValuesByExtrapolation(helpZPreconBounLeft,
                                                                                   helpZPreconBounRight);
                ZPreConBoundaryLeft = PreConBoundary.LeftValue;
                ZPreConBoundaryRight = PreConBoundary.RightValue;
            }
            else
            {
                ZPreConBoundaryLeft = -double.NaN;
                ZPreConBoundaryRight = -double.NaN;
                helpZPreconBounLeft = -double.NaN;
                helpZPreconBounRight = -double.NaN;
            }

            // determine loads if any. A distinction is made between temporary and permanent loads because temporary loads will affect the pore pressures 
            // (by means of the degree of consolidation on loads) but the permanent loads not.
            if (!(uniformLoads == null))
            {
                double xMid = 0.5 * (XLeft + XRight);
                TopLoadPermLeft = 0;
                TopLoadPermRight = 0;
                TopLoadTempLeft = 0;
                TopLoadTempRight = 0;

                foreach (var uniformLoad in uniformLoads)
                {
                    if (uniformLoad.LoadType == UniformLoad.LoadTypeEnum.Permanent)
                    {
                        if ((uniformLoad.XStart < xMid) && (uniformLoad.XEnd > xMid))
                        {
                            TopLoadPermLeft = TopLoadPermLeft + uniformLoad.Magnitude;
                            TopLoadPermRight = TopLoadPermRight + uniformLoad.Magnitude;
                        }
                    }
                    else
                    {
                        if ((uniformLoad.XStart < xMid) && (uniformLoad.XEnd > xMid))
                        {
                            TopLoadTempLeft = TopLoadTempLeft + uniformLoad.Magnitude;
                            TopLoadTempRight = TopLoadTempRight + uniformLoad.Magnitude;
                        }
                    }
                }
            }
        }

        /// <summary>
        /// create interpolation areas by extrapolating boundary values, regarding x values
        /// on the basis of soilprofile2d
        /// </summary>
        /// <param name="profile2D"></param>
        public void CreateAreas(SoilProfile2D profile2D, Waternet waternet,
                                GeometryPointString preloadSurfaceLine, DegreeofConsolidationMatrix consolidationMatrix)
        {
            List<Waternet.IntersectWaterLine> WaterLineListLeft;
            List<Waternet.IntersectWaterLine> WaterLineListRight;

            // find in geometry layers with thickness 
            // boundaries are freatic lines and layer boundaries
            // vertical13.Freaticlevel = MStabDatafunctions.ZTopAtSurface(x13, FreaticLine);

            // waternet lines
            if (waternet != null && waternet.PhreaticLine != null)
            {
                WaterLineListLeft = waternet.IntersectListOnlyWaterlines(helpLeftX);
                WaterLineListRight = waternet.IntersectListOnlyWaterlines(helpRightX);
            }
            else
            {
                WaterLineListLeft = null;
                WaterLineListRight = null;
            }

            if (preloadSurfaceLine != null)
            {
                helpZPreconBounLeft = preloadSurfaceLine.GetZAtX((helpLeftX));
                helpZPreconBounRight = preloadSurfaceLine.GetZAtX((helpRightX));
            }

            helpLeftpoints = getZcoordinates(helpProfileLeft, helpZFreaLeft, WaterLineListLeft, helpZPreconBounLeft);
            helpRightpoints = getZcoordinates(helpProfileRight, helpZFreaRight, WaterLineListRight, helpZPreconBounRight);

            // force the help points to have equal number of points
            // due to rounding they may be different
            if (helpLeftpoints.Count != helpRightpoints.Count)
            {
                List<double> biggestList = helpLeftpoints.Count > helpRightpoints.Count ? helpLeftpoints : helpRightpoints;
                List<double> smallestList = helpLeftpoints.Count < helpRightpoints.Count ? helpLeftpoints : helpRightpoints;

                while (biggestList.Count > smallestList.Count)
                {
                    double smallestDifference = double.MaxValue;
                    int removeIndex = -1;
                    for (int i = 0; i < biggestList.Count - 1; i++)
                    {
                        if (Math.Abs(biggestList[i + 1] - biggestList[i]) < smallestDifference)
                        {
                            removeIndex = i + 1;
                            smallestDifference = Math.Abs(biggestList[i + 1] - biggestList[i]);
                        }
                    }

                    // Only apply when this is really a small difference
                    if (smallestDifference < Constants.CGeoAccu*3)
                    {
                        biggestList.RemoveAt(removeIndex);
                    }
                    else
                    {
                        break;
                    }
                }
            }

            // check if both sides have the same number of layers (areas)
            if (helpLeftpoints.Count != helpRightpoints.Count)
            {
                var translate = this.Translate("NumberOfAreasAtLeftAndRightSidesOfColumnIsDifferent{0}{1}");
                var message = String.Format(translate, XLeft, XRight);
                throw new Exception(message);
            }
            else
            {
                for (int i = 0; i < helpLeftpoints.Count - 1; i++)
                {
                    if (helpLeftpoints[i] <= helpProfileLeft.Layers[0].TopLevel + Constants.CAlmostZero)
                    {
                        var area = new InterpolationArea();
                        area.HelpZTopLeft = helpLeftpoints[i];
                        area.HelpZTopRight = helpRightpoints[i];
                        area.HelpZBottomLeft = helpLeftpoints[i + 1];
                        area.HelpZBottomRight = helpRightpoints[i + 1];

                        if (i == helpLeftpoints.Count - 2) // extra area
                        {
                            area.Soil = areas[areas.Count - 1].Soil; // soil from former area
                            area.SoilLayer = null;
                        }
                        else
                        {
                            area.Soil = getSoilAtZ(helpProfileLeft, 0.5*(area.HelpZTopLeft + area.HelpZBottomLeft));
                            var aPoint = new GeometryPoint(helpLeftX, 0, 0.5*(area.HelpZTopLeft + area.HelpZBottomLeft));

                            area.SoilLayer = profile2D.GetSoilLayer(aPoint);
                        }

                        object loadConsolidator = null;
                        foreach (var consolidator in consolidationMatrix.Consolidators)
                        {
                            if (consolidator is UniformLoad)
                            {
                                loadConsolidator = consolidator;
                                break;
                            }
                        }

                        area.DegreeOfConsolidationDueToLoad = 0;
                        if (loadConsolidator != null && area.SoilLayer != null)
                        {
                            area.DegreeOfConsolidationDueToLoad = consolidationMatrix.GetConsolidationValue(loadConsolidator, area.SoilLayer).Value;
                        }

                        // initialize Pop, phi and cohesion from soil
                        if (area.Soil.UsePop)
                        {
                            area.POP = area.Soil.POP;
                        }
                        else
                        {
                            /*area. =
                            area.Soil.PreConsolidationStressValue;*/
                        }
                        area.Phi = area.Soil.FrictionAngle;
                        area.Cohesion = area.Soil.Cohesion;

                        //TODO: Han check
                        area.ShearStrenghModel = area.Soil.ShearStrengthModel;

                        AreaBoundary topBoundary = GetAreaBoundaryValuesByExtrapolation(area.HelpZTopLeft,
                                                                                        area.HelpZTopRight);
                        AreaBoundary bottomBoundary = GetAreaBoundaryValuesByExtrapolation(area.HelpZBottomLeft,
                                                                                           area.HelpZBottomRight);

                        area.TopLeftPoint = new InterpolationPoint()
                            {
                                XCoordinate = XLeft, ZCoordinate = topBoundary.LeftValue
                            };
                        area.TopRightPoint = new InterpolationPoint()
                            {
                                XCoordinate = XRight, ZCoordinate = topBoundary.RightValue
                            };
                        area.BottomLeftPoint = new InterpolationPoint()
                            {
                                XCoordinate = XLeft, ZCoordinate = bottomBoundary.LeftValue
                            };
                        area.BottomRightPoint = new InterpolationPoint()
                            {
                                XCoordinate = XRight, ZCoordinate = bottomBoundary.RightValue
                            };
                        areas.Add(area);
                    }
                }
            }
        }

        /// <summary>
        /// get extrapolated values, regarding x, of area boundary
        /// </summary>
        /// <param name="valueOneThird"></param>
        /// <param name="valueTwoThird"></param>
        /// <returns></returns>
        internal AreaBoundary GetAreaBoundaryValuesByExtrapolation(double valueOneThird, double valueTwoThird)
        {
            double valueLeft;
            double valueRight;
            double aTanAngle = 0;
            if (Math.Abs(valueOneThird - valueTwoThird) < 0.000001)
            {
                valueLeft = valueOneThird;
                valueRight = valueTwoThird;
                return new AreaBoundary
                    {
                        LeftValue = valueLeft,
                        RightValue = valueRight,
                        AtanAngle = aTanAngle
                    };
            }
            aTanAngle = Math.Atan2((valueTwoThird - valueOneThird), (helpRightX - helpLeftX));
            valueLeft = valueOneThird - Math.Tan(aTanAngle)*(helpLeftX - XLeft);
            valueRight = valueTwoThird + Math.Tan(aTanAngle)*(XRight - helpRightX);
            return new AreaBoundary
                {
                    LeftValue = valueLeft,
                    RightValue = valueRight,
                    AtanAngle = aTanAngle
                };
        }

        internal void SetHelpLeftX(int leftX)
        {
            helpLeftX = leftX;
        }

        internal void SetHelpRightX(int rightX)
        {
            helpRightX = rightX;
        }

        /// <summary>
        /// Check if aZCoor already excists in result
        /// </summary>
        /// <param name="result"></param>
        /// <param name="aZCoor"></param>
        /// <returns></returns>
        private bool isNotPresent(List<double> result, double aZCoor)
        {
            bool notpresent = true;

            foreach (var z in result)
            {
                if (Math.Abs(aZCoor - z) <= Constants.CGeoAccu)
                {
                    notpresent = false;
                    break;
                }
            }
            return notpresent;
        }

        /// <summary>
        /// The area boundaries are on layer separations
        /// On the freatic line (if below surface)
        /// on waternet lines
        /// </summary>
        /// <param name="Profile1d"></param>
        /// <param name="aZFrea"></param>
        /// <param name="aWaterLineList"></param>
        /// <returns></returns>
        private List<double> getZcoordinates(SoilProfile1D Profile1d, double aZFrea, List<Waternet.IntersectWaterLine> aWaterLineList, double helpZPreconBoun)
        {
            // add soil data
            var result = new List<double>();
            for (int i = 0; i < Profile1d.LayerCount; i++)
            {
                if (isNotPresent(result, Profile1d.Layers[i].TopLevel))
                {
                    result.Add(Profile1d.Layers[i].TopLevel);
                }
            }

            if (isNotPresent(result, Profile1d.Layers[Profile1d.LayerCount - 1].BottomLevel))
            {
                result.Add(Profile1d.Layers[Profile1d.LayerCount - 1].BottomLevel);
            }

            if (!double.IsNaN(helpZPreconBoun))
            {
                if (isNotPresent(result, helpZPreconBoun))
                {
                    result.Add(helpZPreconBoun);
                }
            }

            // add freatic levels
            if (isNotPresent(result, aZFrea))
            {
                // ZFrea is only a layer seperation if it is below the surface
                if (aZFrea < Profile1d.Layers[0].TopLevel)
                {
                    result.Add(aZFrea);
                }
            }

            // add water net data
            if (aWaterLineList != null && aWaterLineList.Count > 0)
            {
                for (int i = 0; i < aWaterLineList.Count; i++)
                {
                    if (isNotPresent(result, aWaterLineList[i].Z))
                    {
                        result.Add(aWaterLineList[i].Z);
                    }
                }
            }

            result.Add(Profile1d.Layers[Profile1d.LayerCount - 1].BottomLevel - Constants.CGeometryExtension);

            // sort the z values this will be from top to bottom
            result.Sort();
            result.Reverse();

            return result;
        }

        /// <summary>
        /// returns the soil from a profile at a certain z
        /// </summary>
        /// <param name="aSoilProfile"></param>
        /// <param name="aZCoor"></param>
        /// <returns></returns>
        private Soil getSoilAtZ(SoilProfile1D aSoilProfile, double aZCoor)
        {
            for (int i = 0; i < aSoilProfile.LayerCount; i++)
            {
                if ((aZCoor <= aSoilProfile.Layers[i].TopLevel) && (aZCoor >= aSoilProfile.Layers[i].BottomLevel))
                {
                    return aSoilProfile.Layers[i].Soil;
                }
            }
            return aSoilProfile.Layers[aSoilProfile.LayerCount - 1].Soil;
        }

        /// <summary>
        /// This procedure calculates the soil stress and the total stress (= soil stress + free water above surface line + loads) at the 4 corners of an area.
        /// </summary>
        private void FillPointsWithTotalStressData()
        {
            double deltaStressLeft;
            double deltaStressRight;
            double LSoilStressLeft = 0.0;
            double LSoilStressRight = 0.0;
            double LTotalStressLeft = 0.0;
            double LTotalStressRight = 0.0;
            double waterOnSurfaceAtLeft = 0.0;
            double waterOnSurfaceAtRight = 0.0;
            AreaBoundary soilStress;
            AreaBoundary totalStress;
            AreaBoundary waterOnSurface;
            
            // At the top of the toppest area : the soil stress is zero and the total stress is the free water above surface line
            if (Areas.Count > 0)
            {
                if (helpZFreaLeft > areas[0].HelpZTopLeft)
                {
                    waterOnSurfaceAtLeft = (helpZFreaLeft - areas[0].HelpZTopLeft) * GammaWater;
                }
                if (helpZFreaRight > areas[0].HelpZTopRight)
                {
                    waterOnSurfaceAtRight = (helpZFreaRight - areas[0].HelpZTopRight) * GammaWater;
                }
                waterOnSurface = GetAreaBoundaryValuesByExtrapolation(waterOnSurfaceAtLeft, waterOnSurfaceAtRight);
                
                areas[0].TopLeftPoint.SoilStress = 0.0;
                areas[0].TopRightPoint.SoilStress = 0.0;

                areas[0].TopLeftPoint.TotalStress = waterOnSurface.LeftValue + UniformLoadPermMagnitude + UniformLoadTempMagnitude;
                areas[0].TopRightPoint.TotalStress = waterOnSurface.RightValue + UniformLoadPermMagnitude + UniformLoadTempMagnitude;
            }

            for (int i = 0; i < Areas.Count; i++)
            {
                areas[i].Phi = areas[i].Soil.FrictionAngle;
                areas[i].Cohesion = areas[i].Soil.Cohesion;
                areas[i].DilatancyType = areas[i].Soil.DilatancyType;
                areas[i].POP = areas[i].Soil.POP;
                areas[i].OCR = areas[i].Soil.OCR;

                switch (areas[i].DilatancyType)
                {
                    case DilatancyType.Phi:
                        areas[i].Dilatancy = areas[i].Soil.FrictionAngle;
                        break;
                    case DilatancyType.Zero:
                        areas[i].Dilatancy = 0.0;
                        break;
                    case DilatancyType.MinusPhi:
                        areas[i].Dilatancy = -areas[i].Soil.FrictionAngle;
                        break;
                }
                double LSinPhi = Math.Sin(areas[i].Phi*CDegreeToRad);
                double LCosPhi = Math.Cos(areas[i].Phi*CDegreeToRad);
                double LCosDilatancy = Math.Cos(areas[i].Dilatancy*CDegreeToRad);
                double lSinDilatancy = Math.Sin(areas[i].Dilatancy*CDegreeToRad);

                double LSigTerm = LCosDilatancy*LSinPhi/(1 - lSinDilatancy*LSinPhi);
                double LCohTerm = LCosDilatancy*LCosPhi/(1 - lSinDilatancy*LSinPhi);

                areas[i].Phi = Math.Atan(LSigTerm)*CRadToDegree;
                areas[i].Cohesion = areas[i].Cohesion*LCohTerm;

                // TODO dilatancy must be taken as soil parameter first
                //  case DilatancyType.diUserdefined:
                //  areas[i].Dilatancy = areas[i].Soil.Dilatancy;
                //  break;
                if (i > 0)
                {
                    areas[i].TopLeftPoint.SoilStress = areas[i - 1].BottomLeftPoint.SoilStress;
                    areas[i].TopRightPoint.SoilStress = areas[i - 1].BottomRightPoint.SoilStress;
                    areas[i].TopLeftPoint.TotalStress = areas[i - 1].BottomLeftPoint.TotalStress;
                    areas[i].TopRightPoint.TotalStress = areas[i - 1].BottomRightPoint.TotalStress;
                }

                // calculate total stress
                if (0.5*(areas[i].HelpZTopLeft + areas[i].HelpZBottomLeft) < helpZFreaLeft)
                {
                    // if phreatic line above soil
                    deltaStressLeft = (areas[i].HelpZTopLeft - areas[i].HelpZBottomLeft)*
                                      areas[i].Soil.BelowPhreaticLevel;
                    deltaStressRight = (areas[i].HelpZTopRight - areas[i].HelpZBottomRight)*
                                       areas[i].Soil.BelowPhreaticLevel;
                }
                else
                {
                    deltaStressLeft = (areas[i].HelpZTopLeft - areas[i].HelpZBottomLeft) * areas[i].Soil.AbovePhreaticLevel;
                    deltaStressRight = (areas[i].HelpZTopRight - areas[i].HelpZBottomRight) * areas[i].Soil.AbovePhreaticLevel;
                }
                LSoilStressLeft = LSoilStressLeft + deltaStressLeft;
                LSoilStressRight = LSoilStressRight + deltaStressRight;
                soilStress = GetAreaBoundaryValuesByExtrapolation(LSoilStressLeft, LSoilStressRight);
                areas[i].BottomLeftPoint.SoilStress = soilStress.LeftValue;
                areas[i].BottomRightPoint.SoilStress = soilStress.RightValue;

                LTotalStressLeft = LSoilStressLeft + waterOnSurfaceAtLeft + UniformLoadTempMagnitude + UniformLoadPermMagnitude;
                LTotalStressRight = LSoilStressRight + waterOnSurfaceAtRight + UniformLoadTempMagnitude + UniformLoadPermMagnitude;
                totalStress = GetAreaBoundaryValuesByExtrapolation(LTotalStressLeft, LTotalStressRight);
                areas[i].BottomLeftPoint.TotalStress = totalStress.LeftValue;
                areas[i].BottomRightPoint.TotalStress = totalStress.RightValue;
            }

            // Su measuresd
            int j = 0;

            // the bottom area is added, so don't interpolate over that area
            // give a special treatment
            while (j < Areas.Count - 1)
            {
                if (areas[j].Soil.ShearStrengthModel == ShearStrengthModel.CuMeasured)
                {
                    int startIndex = j;
                    int endIndex = j;

                    while ((endIndex < Areas.Count - 2) && (areas[endIndex + 1].Soil == (areas[startIndex].Soil)))
                    {
                        // there might be more areas with the same soil
                        endIndex++;
                    }

                    double leftZTopCu = areas[startIndex].TopLeftPoint.ZCoordinate;
                    double leftZBottomCu = areas[endIndex].BottomLeftPoint.ZCoordinate;
                    double cuTop = areas[startIndex].Soil.CuTop;
                    double cuBotom = areas[startIndex].Soil.CuBottom;
                    double rightZTopCu = areas[startIndex].TopRightPoint.ZCoordinate;
                    double rightZBottomCu = areas[endIndex].BottomRightPoint.ZCoordinate;

                    for (int k = startIndex; k < endIndex + 1; k++)
                    {
                        // left side
                        areas[k].TopLeftPoint.SuMeasured = MStabDatafunctions.LinInpolY(leftZTopCu, cuTop, leftZBottomCu, cuBotom, areas[k].TopLeftPoint.ZCoordinate);
                        areas[k].BottomLeftPoint.SuMeasured = MStabDatafunctions.LinInpolY(leftZTopCu, cuTop, leftZBottomCu, cuBotom, areas[k].BottomLeftPoint.ZCoordinate);

                        // right side
                        areas[k].TopRightPoint.SuMeasured = MStabDatafunctions.LinInpolY(rightZTopCu, cuTop, rightZBottomCu, cuBotom, areas[k].TopRightPoint.ZCoordinate);
                        areas[k].BottomRightPoint.SuMeasured = MStabDatafunctions.LinInpolY(rightZTopCu, cuTop, rightZBottomCu, cuBotom, areas[k].BottomRightPoint.ZCoordinate);

                        areas[k].TopLayerArea = areas[startIndex];
                        areas[k].BottomLayerArea = areas[endIndex];
                    }

                    j = endIndex;
                }
                j++;
            }
            if (Areas.Count == 0)
            {
                return;
            }
            if (areas[Areas.Count - 1].Soil.ShearStrengthModel == ShearStrengthModel.CuMeasured)
            {
                // left side
                double ZTop = areas[Areas.Count - 2].TopLeftPoint.ZCoordinate;
                double ZBot = areas[Areas.Count - 2].BottomLeftPoint.ZCoordinate;
                double LRatio = 0;
                if ((ZTop - ZBot) > 0)
                {
                    LRatio = (areas[Areas.Count - 2].TopLeftPoint.SuMeasured -
                              areas[Areas.Count - 2].BottomLeftPoint.SuMeasured)/
                             (ZTop - ZBot);
                }
                areas[Areas.Count - 1].TopLeftPoint.SuMeasured = areas[Areas.Count - 2].BottomLeftPoint.SuMeasured;
                areas[Areas.Count - 1].BottomLeftPoint.SuMeasured = areas[Areas.Count - 1].TopLeftPoint.SuMeasured + LRatio*(areas[Areas.Count - 1].TopLeftPoint.ZCoordinate - areas[Areas.Count - 1].BottomLeftPoint.ZCoordinate);

                // Right side
                // left side
                ZTop = areas[Areas.Count - 2].TopRightPoint.ZCoordinate;
                ZBot = areas[Areas.Count - 2].BottomRightPoint.ZCoordinate;
                LRatio = 0;
                if ((ZTop - ZBot) > 0)
                {
                    LRatio = (areas[Areas.Count - 2].TopRightPoint.SuMeasured -
                              areas[Areas.Count - 2].BottomRightPoint.SuMeasured)/
                             (ZTop - ZBot);
                }
                areas[Areas.Count - 1].TopRightPoint.SuMeasured = areas[Areas.Count - 2].BottomRightPoint.SuMeasured;
                areas[Areas.Count - 1].BottomRightPoint.SuMeasured =
                    areas[Areas.Count - 1].TopRightPoint.SuMeasured +
                    LRatio*(areas[Areas.Count - 1].TopRightPoint.ZCoordinate - areas[Areas.Count - 1].BottomRightPoint.ZCoordinate);
            }
        }

        /// <summary>
        /// The piezometric pressures (from PL-lines) are calculated at the 4 corners of an area.
        /// </summary>
        /// <param name="waternet">Waternet</param>
        private void FillPointsWithWaternetPressures(Waternet waternet)
        {
            for (int i = 0; i < Areas.Count; i++)
            {
                double helpPorePrestopLeft = waternet.GetPorePressure(helpLeftX, areas[i].HelpZTopLeft, double.MaxValue, true);
                double helpPorePrestopRight = waternet.GetPorePressure(helpRightX, areas[i].HelpZTopRight, double.MaxValue, true);
                AreaBoundary toppore = GetAreaBoundaryValuesByExtrapolation(helpPorePrestopLeft, helpPorePrestopRight);
                areas[i].TopLeftPoint.PnLinePore = toppore.LeftValue;
                areas[i].TopRightPoint.PnLinePore = toppore.RightValue;

                double helpPorePresbotLeft = waternet.GetPorePressure(helpLeftX, areas[i].HelpZBottomLeft, double.MaxValue, false);
                double helpPorePresbotRight = waternet.GetPorePressure(helpRightX, areas[i].HelpZBottomRight, double.MaxValue, false);
                AreaBoundary botpore = GetAreaBoundaryValuesByExtrapolation(helpPorePresbotLeft, helpPorePresbotRight);
                areas[i].BottomLeftPoint.PnLinePore = botpore.LeftValue;
                areas[i].BottomRightPoint.PnLinePore = botpore.RightValue;
            }
        }

        /// <summary>
        /// The extra pore pressure caused by temporary loads are calculated at the 4 corners of an area.
        /// </summary>
        private void FillPointsWithPorePressuresDueToDegreeOfConsolidationLoad()
        {
            for (int i = 0; i < Areas.Count; i++)
            {
                // if above phreatic line, temporary load is the same as permanent load (no excess pore pressure)
                double tempLoadStress;
                double zMiddleLeft = 0.5 * (areas[i].TopLeftPoint.ZCoordinate + areas[i].BottomLeftPoint.ZCoordinate);
                double zMiddleRight = 0.5 * (areas[i].TopRightPoint.ZCoordinate + areas[i].BottomRightPoint.ZCoordinate);
                if ((zMiddleLeft > ZFreaLeft) && (zMiddleRight > ZFreaRight))
                {
                    tempLoadStress = 0;
                }
                else
                {
                    tempLoadStress = UniformLoadTempMagnitude;
                }
                double porePressure = (1 - areas[i].DegreeOfConsolidationDueToLoad/100.0) * tempLoadStress;
                areas[i].TopLeftPoint.PorePressureDueToDegreeofConsolidationLoad = porePressure;
                areas[i].TopRightPoint.PorePressureDueToDegreeofConsolidationLoad = porePressure;
                areas[i].BottomLeftPoint.PorePressureDueToDegreeofConsolidationLoad = porePressure;
                areas[i].BottomRightPoint.PorePressureDueToDegreeofConsolidationLoad = porePressure;
            }
        }

        private void FillPointsWithTotalPorePressures()
        {
            for (int i = 0; i < Areas.Count; i++)
            {
                areas[i].TopLeftPoint.TotalPore = areas[i].TopLeftPoint.PnLinePore + areas[i].TopLeftPoint.PorePressureDueToDegreeofConsolidationLoad;
                areas[i].TopRightPoint.TotalPore = areas[i].TopRightPoint.PnLinePore + areas[i].TopRightPoint.PorePressureDueToDegreeofConsolidationLoad;
                areas[i].BottomLeftPoint.TotalPore = areas[i].BottomLeftPoint.PnLinePore + areas[i].BottomLeftPoint.PorePressureDueToDegreeofConsolidationLoad;
                areas[i].BottomRightPoint.TotalPore = areas[i].BottomRightPoint.PnLinePore + areas[i].BottomRightPoint.PorePressureDueToDegreeofConsolidationLoad;
            }
        }

        private void FillPointsWithEffectiveStresses()
        {
            for (int i = 0; i < Areas.Count; i++)
            {
                areas[i].TopLeftPoint.EffectiveStress = areas[i].TopLeftPoint.TotalStress - areas[i].TopLeftPoint.TotalPore;
                areas[i].TopRightPoint.EffectiveStress = areas[i].TopRightPoint.TotalStress - areas[i].TopRightPoint.TotalPore;
                areas[i].BottomLeftPoint.EffectiveStress = areas[i].BottomLeftPoint.TotalStress - areas[i].BottomLeftPoint.TotalPore;
                areas[i].BottomRightPoint.EffectiveStress = areas[i].BottomRightPoint.TotalStress - areas[i].BottomRightPoint.TotalPore;
            }
        }

        private void FillPointsWithPreLoadStressData()
        {
            double deltaPreLoadStressLeft;
            double deltaPreLoadStressRight;
            double LTotPreloadStressLeft = 0.0;
            double LTotPreloadStressRight = 0.0;

            if (Areas.Count > 0)
            {
                areas[0].TopLeftPoint.PreLoadStress = 0.0;
                areas[0].TopRightPoint.PreLoadStress = 0.0;
            }

            for (int i = 0; i < Areas.Count; i++)
            {
                if (i > 0)
                {
                    areas[i].TopLeftPoint.PreLoadStress = areas[i - 1].BottomLeftPoint.PreLoadStress;
                    areas[i].TopRightPoint.PreLoadStress = areas[i - 1].BottomRightPoint.PreLoadStress;
                }

                // calulate preloadstress
                if (ZPreConBoundaryLeft != double.NaN)
                {
                    if (areas[i].HelpZTopLeft <= ZPreConBoundaryLeft + Constants.CGeoAccu)
                    {
                        if (0.5*(areas[i].HelpZTopLeft + areas[i].HelpZBottomLeft) < helpZFreaLeft)
                        {
                            deltaPreLoadStressLeft = (areas[i].HelpZTopLeft - areas[i].HelpZBottomLeft)*
                                                     areas[i].Soil.BelowPhreaticLevel;
                            deltaPreLoadStressRight = (areas[i].HelpZTopRight - areas[i].HelpZBottomRight)*
                                                      areas[i].Soil.BelowPhreaticLevel;
                        }
                        else
                        {
                            deltaPreLoadStressLeft = (areas[i].HelpZTopLeft - areas[i].HelpZBottomLeft)*
                                                     areas[i].Soil.AbovePhreaticLevel;
                            deltaPreLoadStressRight = (areas[i].HelpZTopRight - areas[i].HelpZBottomRight)*
                                                      areas[i].Soil.AbovePhreaticLevel;
                        }
                        LTotPreloadStressLeft = LTotPreloadStressLeft + deltaPreLoadStressLeft;
                        LTotPreloadStressRight = LTotPreloadStressRight + deltaPreLoadStressRight;

                        AreaBoundary PreLoadStress = GetAreaBoundaryValuesByExtrapolation(LTotPreloadStressLeft,
                                                                                          LTotPreloadStressRight);

                        areas[i].BottomLeftPoint.PreLoadStress = PreLoadStress.LeftValue;
                        areas[i].BottomRightPoint.PreLoadStress = PreLoadStress.RightValue;
                    }
                    else
                    {
                        areas[i].BottomLeftPoint.PreLoadStress = 0;
                        areas[i].BottomRightPoint.PreLoadStress = 0;
                    }
                }
            }
        }

        /// <summary>
        /// Round double to 3 digits
        /// </summary>
        /// <param name="valueLeft"></param>
        /// <returns></returns>
        private static double RoundThreeDigits(double valueLeft)
        {
            return Math.Round(valueLeft, 3, MidpointRounding.AwayFromZero);
        }

        /// <summary>
        /// get x value at 2/3 of the total value
        /// </summary>
        /// <returns></returns>
        private double GetXAtTwoThird()
        {
            var result = XLeft + 2*(XRight - XLeft)/3;
            return result;
        }

        /// <summary>
        /// get x value at 1/3 of the total value
        /// </summary>
        /// <returns></returns>
        private double GetXAtOneThird()
        {
            var result = XLeft + (XRight - XLeft)/3;
            return result;
        }

        public override string ToString()
        {
            return string.Format("{0:F2} - {1:F2}", XLeft, XRight);
        }
    }
}