using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Consolidation;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Standard;
using Deltares.Standard.Language;
using Deltares.Standard.Logging;

// =======================================================================================================================
// Class name:PreProcesStabilityCalculation
//
// Description:  Executes the calculation of slice values
// at fixed points in geometry
//
// If this data is needed in the slices it can be accesed by interpolation
//  
// First all the points in the geometry (layer points and pl-line points) 
// waternet line points, intersection points from waternetlines with layers
// and intersection point freatic line with surface and layers are searched for
// the start and end point of a traffic load is a point as well
// After sorting columns are made going trough al these points
// Each column is subdivided in areas  (soil layers) where the freatic line is a boundary as well
// with four points;
// The lineair data is calculated in those points
//
//            column1         colum2
//       |-------------------------------\
//       |p1           p2 |p1         p2 |\  column 3
//       |                |              |p1\      
//       |                |              |    \    
//       |                |              |     \   
//       |    area 0      |   area 0     |    p2 \    as the freatic line is a boundary 
//       |                |              |       |    the soil in an  area is either dry or wet
//       |                |              | area0 |    
//       |                |              |       |    top of area 0 is the surface 
//       |p3            p4| p3         p4|p3  p4 |    for each column there is column data
//       ----------------------------------------|                                        
//       |p1            p2|p1          p2|p1   p2|    angle on surface
//       |                |    x z1      |       |    (water) load on surface 
//       |                |              |       |    this data is needed for each slice
//       |                |      area 1  | area 1|  
//       |    area 1      |              |       |    all other data is calculated in the points p1,p2 etc
//       |                |              |       |    for each area
//       |p3            p4|p3          p4|p3   p4|  
//       ----------------------------------------|    to find data during calculation for point z1 then                                     
//                                                    first the column is found (column 2)
//                                                    next the area (area 1) and data is found 
//                                                    by lineair interpolation between the points p1 t/m p4
//
// Copyright (c) 2008-2013 Deltares
//
// Date        ID   Modification
// 2013-03-27  Best Created
// ======================================================================================================================

namespace Deltares.Stability.Calculation2
{
    public class PreProcesStabilityCalculation
    {
        internal List<CalculatePopSoil> popsInSoilArray = new List<CalculatePopSoil>();
        internal List<double> xCoordList = new List<double>();
        private List<GeometryPoint> changingPoints = new List<GeometryPoint>();
        private SlipplaneConstraints constraints;
        private DegreeofConsolidationMatrix degreeOfConsolidationLoads = new DegreeofConsolidationMatrix();
        private List<InterpolationColumn> interpolationColumns = new List<InterpolationColumn>();
        private int lastUsedColumn = 0;

        //  ... Data needed for calculation of pop from preconsolidation stress

        private List<MultiplicationFactorOnCPhiForUplift> multiplicationFactorsArray = new List<MultiplicationFactorOnCPhiForUplift>();
        private List<PreConsolidationStress> preconStressArray = new List<PreConsolidationStress>();
        private double preconsolidationStressModelFactor = 1;
        private GeometryPointString preloadSurfaceLine = new GeometryPointString();
        private SoilProfile2D profile2D;
        private SurfaceLine2 surfaceLine;
        private List<UniformLoad> uniformLoads = new List<UniformLoad>();
        private Waternet waternet;

        // ... Only for debugging 
//        public StreamWriter writer = new StreamWriter("test.txt");


        //  ... End of data needed for calculation of pop from preconsolidation stress

        public SlipplaneConstraints Constraints
        {
            get
            {
                return constraints;
            }
            //set { constraints = value; }
        }

        public SoilProfile2D Profile2D
        {
            get
            {
                return profile2D;
            }
        }

        public Waternet Waternet
        {
            get
            {
                return waternet;
            }
        }

        public SurfaceLine2 SurfaceLine
        {
            get
            {
                return surfaceLine;
            }
        }

        public List<UniformLoad> UniformLoads
        {
            get
            {
                return uniformLoads;
            }
        }

        public List<InterpolationColumn> InterpolationColumns
        {
            get
            {
                return interpolationColumns;
            }
        }

        public List<GeometryPoint> ChangingPoints
        {
            get
            {
                return changingPoints;
            }
        }

        public GeometryPointString PreloadSurfaceLine
        {
            get
            {
                return preloadSurfaceLine;
            }
        }

        /// <summary>
        /// Model factor for preconsolidation stress
        /// </summary>
        public double PreconsolidationStressModelFactor
        {
            get
            {
                return preconsolidationStressModelFactor;
            }
            set
            {
                preconsolidationStressModelFactor = value;
            }
        }

        public void SetSoilProfile2D(SoilProfile2D profile)
        {
            profile2D = profile;
        }

        public void SetPreConsolidationStress(List<PreConsolidationStress> preConsolidationStressList)
        {
            preconStressArray = preConsolidationStressList;
        }

        public void SetMultiplicationFactorOnCPhi(List<MultiplicationFactorOnCPhiForUplift> multiplicationFactorOnCPhiList)
        {
            multiplicationFactorsArray = multiplicationFactorOnCPhiList;
        }

        public void SetWaternet(Waternet waternetExt)
        {
            waternet = waternetExt;
        }

        public void SetSurfaceLine(SurfaceLine2 surfaceLine)
        {
            this.surfaceLine = surfaceLine;
        }

        public void SetUniformLoads(List<UniformLoad> uniformLoads)
        {
            this.uniformLoads = uniformLoads;
        }

        public void SetDegreeOfConsolidationLoads(DegreeofConsolidationMatrix consolidationMatrix)
        {
            degreeOfConsolidationLoads = consolidationMatrix;
        }

        public void SetPreloadSurfaceLine(GeometryPointString preloadSurfaceLine)
        {
            this.preloadSurfaceLine = preloadSurfaceLine;
        }

        /// <summary>
        /// For probabilistic calculation is might be necessary 
        /// to fill the interpolation columns with diffreernt data.
        /// The columns might change due to change in freatic lines 
        /// So the columns has to be newly created
        /// furthermore it is asumed that new unit weights and so on are changed in the soils
        /// </summary>
        public void UpdatePreprocessing()
        {
            interpolationColumns = new List<InterpolationColumn>();
            StartPreProcessing();
        }

        public void SortPoints()
        {
            xCoordList.Sort();
        }

        /// <summary>
        /// Create  a column between all the x coordinates found
        /// </summary>
        public void CreateInterpolationColumns()
        {
            interpolationColumns.Clear();

            for (int i = 0; i < xCoordList.Count - 1; i++)
            {
                var interpolationColumn = new InterpolationColumn();
                interpolationColumn.XLeft = xCoordList[i];
                interpolationColumn.XRight = xCoordList[i + 1];
                interpolationColumns.Add(interpolationColumn);
            }
        }

        /// <summary>
        /// Circels(partly  outside the geometry must be calculates as well. For that reason the column data is extended to both sides of the geomtry
        /// </summary>
        public void AddExtendedInterpolationColumns()
        {
            if (interpolationColumns.Count > 0)
            {
                // extension on the left side of the geometry
                var interpolationColumnL = new InterpolationColumn(interpolationColumns[0]);
                interpolationColumnL.XRight = interpolationColumns[0].XLeft;
                interpolationColumnL.XLeft = interpolationColumnL.XRight - Constants.CGeometryExtension;
                for (int i = 0; i < interpolationColumnL.Areas.Count; i++)
                {
                    interpolationColumnL.Areas[i].POP = interpolationColumns[0].Areas[i].POP;
                    interpolationColumnL.Areas[i].OCR = interpolationColumns[0].Areas[i].OCR;
//                    interpolationColumnL.Areas[i].Cohesion = interpolationColumns[0].Areas[i].Cohesion;
//                    interpolationColumnL.Areas[i].Phi = interpolationColumns[0].Areas[i].Phi;
//                    interpolationColumnL.Areas[i].Dilatancy = interpolationColumns[0].Areas[i].Dilatancy;
//                    interpolationColumnL.Areas[i].Su = interpolationColumns[0].Areas[i].Su;
                    interpolationColumnL.Areas[i].TopRightPoint.CopyStressdata(interpolationColumnL.Areas[i].TopLeftPoint);
                    interpolationColumnL.Areas[i].BottomRightPoint.CopyStressdata(interpolationColumnL.Areas[i].BottomLeftPoint);
                }
                interpolationColumns.Insert(0, interpolationColumnL);

                // extension on the right side of the geometry
                var interpolationColumnR = new InterpolationColumn(interpolationColumns[interpolationColumns.Count - 1]);
                interpolationColumnR.XRight = interpolationColumns[interpolationColumns.Count - 1].XRight + Constants.CGeometryExtension;
                interpolationColumnR.XLeft = interpolationColumnR.XRight - Constants.CGeometryExtension;
                for (int i = 0; i < interpolationColumnR.Areas.Count; i++)
                {
                    interpolationColumnR.Areas[i].POP = interpolationColumns[interpolationColumns.Count - 1].Areas[i].POP;
                    interpolationColumnR.Areas[i].OCR = interpolationColumns[interpolationColumns.Count - 1].Areas[i].OCR;
                    //interpolationColumnL.Areas[i].Cohesion = interpolationColumns[interpolationColumns.Count - 1].Areas[i].Cohesion;
                    //interpolationColumnL.Areas[i].Phi = interpolationColumns[interpolationColumns.Count - 1].Areas[i].Phi;
                    //interpolationColumnL.Areas[i].Dilatancy = interpolationColumns[interpolationColumns.Count - 1].Areas[i].Dilatancy;
                    //interpolationColumnL.Areas[i].Su = interpolationColumns[interpolationColumns.Count - 1].Areas[i].Su;
                    interpolationColumnR.Areas[i].TopLeftPoint.CopyStressdata(interpolationColumnR.Areas[i].TopRightPoint);
                    interpolationColumnR.Areas[i].BottomLeftPoint.CopyStressdata(interpolationColumnR.Areas[i].BottomRightPoint);
                }
                interpolationColumns.Add(interpolationColumnR);
            }
        }

        /// <summary>
        /// As no data is available for the extended columns these are removed first and added later as the data from the other columns is calculated
        /// </summary>
        public void RemoveExtendedInterpolationColumns()
            //===============================================
        {
            if (interpolationColumns.Count > 0)
            {
                interpolationColumns.Remove(interpolationColumns[interpolationColumns.Count - 1]);
                interpolationColumns.Remove(interpolationColumns[0]);
            }
        }

        public void CreateHelpVerticals()
        {
            for (int i = 0; i < interpolationColumns.Count; i++)
            {
                interpolationColumns[i].CreateHelpverticals(profile2D);
            }
        }

        public void AddColumnData()
        {
            for (int i = 0; i < interpolationColumns.Count; i++)
            {
                interpolationColumns[i].AddColumnData(waternet, preloadSurfaceLine, uniformLoads);
            }
        }

        public void CreateAreas()
        {
            for (int i = 0; i < interpolationColumns.Count; i++)
            {
                interpolationColumns[i].CreateAreas(profile2D, waternet, preloadSurfaceLine, degreeOfConsolidationLoads);
            }
        }

        public void DetermineColumnData()
        {
            for (int i = 0; i < interpolationColumns.Count; i++)
            {
                interpolationColumns[i].DetermineColumnData();
            }
        }

        public void FillPointsWithData()
        {
            for (int i = 0; i < interpolationColumns.Count; i++)
            {
                interpolationColumns[i].FillPointsWithData(waternet);
            }
        }

        public InterpolationArea GetAreaAtXZ(double xCoor, double zCoor)
        {
            InterpolationColumn column = GetColumnAtX(xCoor);
            InterpolationArea area = column.getArea(zCoor);
            return area;
        }

        /// <summary>
        /// Get the column with the X in it. Keep track of the last used column to increase search speed.
        /// </summary>
        /// <param name="aXCoor"></param>
        /// <returns></returns>
        public InterpolationColumn GetColumnAtX(double aXCoor)
        {
            if (aXCoor < interpolationColumns[lastUsedColumn].XLeft)
            {
                lastUsedColumn = 0;
            }
            for (int j = lastUsedColumn; j < interpolationColumns.Count; j++)
            {
                if (interpolationColumns[j].XIsInColumn(aXCoor))
                {
                    lastUsedColumn = j;
                    return interpolationColumns[j];
                }
            }
            return null;
        }

        public double DetermineZFreaAtX(double x)
        {
            return GetColumnAtX(x).ZFrea;
        }

        /// <summary>
        /// Determine POP value in an area, only when the soil name of this area is present in the POP list (i.e popsInSoilArray).
        /// If the soil name is present several time within a column, the average is used.
        /// If not, the closest POP value connected to the soil name is used; if several values are found (i.e. several POP values along a vertical), the average 
        /// is used.
        /// </summary>
        /// <param name="xmid">X at middle of column k</param>
        /// <param name="xleft">X at left boundary of column k</param>
        /// <param name="xright">X at left boundary of column k</param>
        /// <param name="popSoils">List of index(es) in the array popsInSoilArray corresponding to the soilName (name may be present more than once). </param>
        /// <returns>POP</returns>
        private double DeterminePop(double xmid, double xleft, double xright, List<CalculatePopSoil> popSoils)
        {
            if (popSoils.Count == 1)
            {
                return popSoils[0].Pop;
            }
            else
            {
                double result = 0;
                double som = 0;
                double number = 0;
                double minDistance = 100000.0;
                bool atLeastOnePopInsideColumn = false;
                for (int i = 0; i < popSoils.Count; i++)
                {
                   if ((popSoils[i].XCoor > xleft) && (popSoils[i].XCoor < xright))
                   {
                        som = som + popSoils[i].Pop;
                        atLeastOnePopInsideColumn = true;
                        number++;
                    }
                    else
                    {
                        double distance = Math.Abs(xmid - popSoils[i].XCoor);
                        if (distance < minDistance)
                        {
                            result = popSoils[i].Pop;
                            minDistance = distance;
                        }
                    }
                }
                if (atLeastOnePopInsideColumn)
                {
                    return som / number;
                }

                // In case the closest value is used, several POP values can be at equi-distance
                for (int i = 0; i < popSoils.Count; i++)
                {
                    double distance = Math.Abs(xmid - popSoils[i].XCoor);
                    if (Math.Abs(distance - minDistance) < 0.001)
                    {
                        som = som + popSoils[i].Pop;
                        number++;
                    }
                }

                if (number > 0)
                {
                    result = som / number;
                }
                
                return result;
            }
        }

        /// <summary>
        /// Fill some areas of column k with POP values. 
        /// Only the areas with a soil name present in the array popsInSoilArray are filled.
        /// </summary>
        /// <param name="k">Column index</param>
        /// <param name="soil">Soil name (from array popsInSoilArray) considered.</param>
        /// <param name="nameIndex">List of index(es) in the array popsInSoilArray corresponding to the soilName (name may be present more than once). </param>
        private void FillAreasInColumnsWithPop(InterpolationColumn interpolationCol, Soil soil, List<CalculatePopSoil> popSoils)
        {
            // shortcuts
            double xLeft = interpolationCol.XLeft;
            double xRight = interpolationCol.XRight;
            double xMid = (xLeft + xRight)*0.5;

            foreach (var area in interpolationCol.Areas)
            {
                if (area.Soil == soil)
                {
                    area.POP = DeterminePop(xMid, xLeft, xRight, popSoils);
                }
            }
        }

        /// <summary>
        /// First determine the list of index(es) in the array popsInSoilArray corresponding to the soilName (name may be present more than once).
        /// Then fill the areas within the column with POP values.
        /// </summary>
        /// <param name="soil">Soil material in the considered area</param>
        public void DeterminePopForNewSoilName(Soil soil)
        {
            List<CalculatePopSoil> popSoils = new List<CalculatePopSoil>(); // name may be present more than once

            for (int k = 0; k < popsInSoilArray.Count; k++)
            {
                if (soil == popsInSoilArray[k].Soil)
                {
                    popSoils.Add(popsInSoilArray[k]);
                }
            }
            // all pre-consolidation stress indexes present in soilName are present in nameIndex
            // now put the POP in all areas with soilName
            // if more than 1 index present, choose the closest to the middle of the column
            // if more than 1 index present in the same column take the average of those
            for (int k = 0; k < interpolationColumns.Count; k++)
            {
                FillAreasInColumnsWithPop(interpolationColumns[k], soil, popSoils);
            }
        }

        public double InterpolatePop(InterpolationColumn Col, int StartIndex, int EndIndex, int searchindex)
            //======================================================================================================
        {
            double zStart = Col.Areas[StartIndex].DetermineZCoorMid();
            double popStart = Col.Areas[StartIndex].POP;
            double zEnd = Col.Areas[EndIndex].DetermineZCoorMid();
            double popEnd = Col.Areas[EndIndex].POP;
            double zSearch = Col.Areas[searchindex].DetermineZCoorMid();
            return MStabDatafunctions.LinInpolY(zStart, popStart, zEnd, popEnd, zSearch);
        }

        /// <summary>
        /// See if a calculated POP is present in this column. If so, calculate the pops of the other areas by interpolation.
        /// If not, leave them for the time beiing. They will be filled in procedure FillNotUsedColumnsWithPop.
        /// </summary>
        /// <param name="k">Column index</param>
        public void FillNotUsedAreasInColumnsWithPopByInterpolation(InterpolationColumn interpolationCol)
        {
            var areaindexWithPop = new List<int>();
            for (int i = 0; i < interpolationCol.Areas.Count; i++)
            {
                if (interpolationCol.Areas[i].POP > -0.001)
                {
                    areaindexWithPop.Add(i);
                }
            }

            if (areaindexWithPop.Count > 0)
            {
                for (int i = 0; i < interpolationCol.Areas.Count; i++)
                {
                    if (interpolationCol.Areas[i].POP < -0.001)
                    {

                        if (i > areaindexWithPop[areaindexWithPop.Count - 1])
                        {
                            interpolationCol.Areas[i].POP = interpolationCol.Areas[areaindexWithPop[areaindexWithPop.Count - 1]].POP;
                        }
                        else
                        {
                            if (i < areaindexWithPop[0])
                            {
                                interpolationCol.Areas[i].POP = interpolationCol.Areas[areaindexWithPop[0]].POP;
                            }
                            else
                            {
                                for (int j = 0; j < areaindexWithPop.Count; j++)
                                {
                                    if ((i > areaindexWithPop[j]) && (i < areaindexWithPop[j + 1]))
                                    {
                                        interpolationCol.Areas[i].POP = InterpolatePop(interpolationCol, areaindexWithPop[j], areaindexWithPop[j + 1], i);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }

        public void StartPreProcessing()
            //===============================
        {
            // temporary data for building purposes
            DeterminePointsAtXFromProfile2DSurfaces();
            SortPoints();

            CreateInterpolationColumns();

            // help verticals are needed to clculate results in the corners of the areas
            CreateHelpVerticals();
            AddColumnData();

            CreateAreas();

            DetermineColumnData();
            FillPointsWithData();

            if (preconStressArray.Count > 0)
            {
                FillAreasWithPreconsolidationStressData();
                // It is possible to use the Preconsolidation Stresses table for some materials (Use POP unchecked) 
                // and to use the given POP for other materials (Use POP checked).
                FillAreasInColumnsWithUserDefinedPop();
            }

            FillAreasInColumnsWithMultiplicationFactor();

            // if Xentry or Xexit lies outside the geometry, then the geometry is extended
            AddExtendedInterpolationColumns();
        }

        public void SetSlipPlaneConstraints(SlipplaneConstraints slipPlaneConstraints)
        {
            constraints = slipPlaneConstraints;
        }

        /// <summary>
        /// Areas are filled in with multiplication factor at left and right. The multiplication factor is applied to reduce the cohesion and the friction angle
        /// (to take uplift into account) for area situated at the right side of the DikeTop At Polder Side and for the layers with a thickness of maximum 4 m
        /// </summary>
        public void FillAreasInColumnsWithMultiplicationFactor()
        {
            GeometryPoint dikeTopAtPolder = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            for (int k = 0; k < interpolationColumns.Count; k++)
            {
                for (int i = 0; i < interpolationColumns[k].Areas.Count; i++)
                {
                    interpolationColumns[k].Areas[i].MultiplicationFactorCPhiLeft = 1.0;
                    interpolationColumns[k].Areas[i].MultiplicationFactorCPhiRight = 1.0;

                    double xLeft = interpolationColumns[k].XLeft;
                    double xRight = interpolationColumns[k].XRight;

                    // If no aquifer present, no reduction applies 
                    if ((profile2D.GetSoilProfile1D(xLeft).GetAquiferLayers().Count > 0) || (profile2D.GetSoilProfile1D(xRight).GetAquiferLayers().Count > 0))
                    {
                        // Determine the thickness of the layer(s) situated above the highest aquifer
                        double thicknessAboveAquiferLeft = profile2D.GetSoilProfile1D(xLeft).TopLevel - profile2D.GetSoilProfile1D(xLeft).GetTopLevelOfHighestAquifer();
                        double thicknessAboveAquiferRight = profile2D.GetSoilProfile1D(xRight).TopLevel - profile2D.GetSoilProfile1D(xRight).GetTopLevelOfHighestAquifer();

                        // Reduction applies only for the points situated at the right side of charac. point DikeTopAtPolderSide
                        if (interpolationColumns[k].XRight > dikeTopAtPolder.X)
                        {
                            // Reduction applies only for a thickness less than 4 m
                            var upliftCalculator = new UpliftCalculator();
                            upliftCalculator.SoilProfile2D = profile2D;
                            upliftCalculator.Waternet = waternet;
                            if ((thicknessAboveAquiferLeft <= 4.0) &&
                                (interpolationColumns[k].Areas[i].BottomLeftPoint.ZCoordinate >= profile2D.GetSoilProfile1D(xLeft).GetTopLevelOfHighestAquifer()))
                            {
                                double upliftFactorLeft = upliftCalculator.DetermineUpliftPotentialAtXAtTopOfHighestAquiferFromA2DProfile(xLeft);
                                interpolationColumns[k].Areas[i].MultiplicationFactorCPhiLeft = GetMultiplicationFactorOnCohAndPhiFromUpliftFactor(upliftFactorLeft);
                            }

                            if ((thicknessAboveAquiferRight <= 4.0) &&
                                (interpolationColumns[k].Areas[i].BottomRightPoint.ZCoordinate >= profile2D.GetSoilProfile1D(xRight).GetTopLevelOfHighestAquifer()))
                            {
                                double upliftFactorRight = upliftCalculator.DetermineUpliftPotentialAtXAtTopOfHighestAquiferFromA2DProfile(xRight);
                                interpolationColumns[k].Areas[i].MultiplicationFactorCPhiRight = GetMultiplicationFactorOnCohAndPhiFromUpliftFactor(upliftFactorRight);
                            }
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Get the multiplication factor (less to 1) which is applied on cohesion and friction angle, depending on the uplift factor, using a
        /// user-defined table of values
        /// </summary>
        /// <param name="upliftFactor">Uplift factor (i.e. mass of soil and free water divided by the water pressures due to PL-lines) at the bottom of the aquitard</param>
        /// <returns></returns>
        public double GetMultiplicationFactorOnCohAndPhiFromUpliftFactor(double upliftFactor)
        {
            double interpolatedMultiplicationFactor = 1.0;
            int number = multiplicationFactorsArray.Count;
            if (multiplicationFactorsArray.Count == 0)
            {
                return 1.0;
            }
            // table is sorted: first value is the minimum Multiplication Factor 
            List<MultiplicationFactorOnCPhiForUplift> sortedMultiplicationFactorsArray = multiplicationFactorsArray.OrderBy(m => m.UpliftFactor).ToList();

            // the uplift factor is more than the maximum uplift factor given by user
            if (upliftFactor >= sortedMultiplicationFactorsArray[number - 1].UpliftFactor)
            {
                return sortedMultiplicationFactorsArray[number - 1].MultiplicationFactor;
            }
            // the uplift factor is less than the minimum uplift factor given by user
            if (upliftFactor <= sortedMultiplicationFactorsArray[0].UpliftFactor)
            {
                return sortedMultiplicationFactorsArray[0].MultiplicationFactor;
            }
            for (int i = 1; i < number; i++)
            {
                if ((upliftFactor < sortedMultiplicationFactorsArray[i].UpliftFactor) &&
                    (upliftFactor > sortedMultiplicationFactorsArray[i - 1].UpliftFactor))
                {
                    double x1 = sortedMultiplicationFactorsArray[i - 1].UpliftFactor;
                    double x2 = sortedMultiplicationFactorsArray[i].UpliftFactor;
                    double y1 = sortedMultiplicationFactorsArray[i - 1].MultiplicationFactor;
                    double y2 = sortedMultiplicationFactorsArray[i].MultiplicationFactor;
                    interpolatedMultiplicationFactor = (y1 - y2)/(x1 - x2)*(upliftFactor - x1) + y1;
                }
            }
            return interpolatedMultiplicationFactor;
        }

        /// <summary>
        /// Coordinates are added to the list where columns are to be made
        /// If a coordinate is les then 1 cm (CVerticalAccu) from an 
        /// excisting vertical, no new point is created
        /// </summary>
        /// <param name="xCoord"></param>
        internal void AddXcoordinate(double xCoord)
        {
            bool present = false;

            foreach (var x in xCoordList)
            {
                if (Math.Abs(x - xCoord) <= Constants.CVerticalAccu)
                {
                    present = true;
                    break;
                }
            }
            if (!present)
            {
                xCoordList.Add(xCoord);
            }
        }

        /// <summary>
        ///  All the changing points in the geometry are searched for
        ///  This includes intersection from freatic line with surface line
        ///  and soil layers
        ///  changing points in piezo lines and waternet lines
        ///  points form layers
        /// </summary>
        internal void DeterminePointsAtXFromProfile2DSurfaces()
        {
            // geometry
            foreach (var surface in Profile2D.Surfaces)
            {
                foreach (var point in surface.GeometrySurface.OuterLoop.Points)
                {
                    AddXcoordinate(point.X);
                    changingPoints.Add(point);
                }
            }
            // waternet

            if (Waternet != null)
            {
                foreach (var waternetline in Waternet.WaternetLineList)
                {
                    foreach (var point in waternetline.Points)
                    {
                        AddXcoordinate(point.X);
                        changingPoints.Add(point);
                    }
                }
                if (Waternet.PhreaticLine != null)
                {
                    foreach (var point in Waternet.PhreaticLine.Points)
                    {
                        AddXcoordinate(point.X);
                        changingPoints.Add(point);
                    }
                }

                foreach (var headline in Waternet.HeadLineList)
                {
                    foreach (var point in headline.Points)
                    {
                        AddXcoordinate(point.X);
                        changingPoints.Add(point);
                    }
                }
            }

            // pre load stress line
            foreach (var point in preloadSurfaceLine.Points)
            {
                AddXcoordinate(point.X);
                changingPoints.Add(point);
            }

            foreach (var uniformLoad in uniformLoads)
            {
                AddXcoordinate(uniformLoad.XStart);
                AddXcoordinate(uniformLoad.XEnd);
                uniformLoad.SurfaceLine2 = SurfaceLine;
            }

            if (uniformLoads.Count == 0)
            {
                // start and end point uniform loads (traffic load in RTO)
                if (SurfaceLine != null &&
                    SurfaceLine.IsDefined(CharacteristicPointType.TrafficLoadInside) &&
                    SurfaceLine.IsDefined(CharacteristicPointType.TrafficLoadOutside))
                {
                    GeometryPoint trafficLoadIn = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.TrafficLoadInside);
                    GeometryPoint trafficLoadOut = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.TrafficLoadOutside);
                    AddXcoordinate(trafficLoadIn.X);
                    AddXcoordinate(trafficLoadOut.X);
                    changingPoints.Add(trafficLoadIn);
                    changingPoints.Add(trafficLoadOut);
                }
            }

            // intersection point(s) freatic line - layers
            // LayerIntersectionPoint
            foreach (var surface in Profile2D.Surfaces)
            {
                if (Waternet != null && Waternet.PhreaticLine != null)
                {
                    IList<GeometryPoint> intersectPoints =
                        surface.GeometrySurface.OuterLoop.ClosedGeometryIntersectionXzPointsWithGeometryPointList(
                            Waternet.PhreaticLine.Points);
                    for (int i = 0; i < intersectPoints.Count; i++)
                    {
                        AddXcoordinate(intersectPoints[i].X);
                        changingPoints.Add(intersectPoints[i]);
                    }
                }
            }

            // intersection point(s) between waternet line(s) (except the surface line) and layers
            foreach (var surface in Profile2D.Surfaces)
            {
                if (Waternet != null)
                {
                    foreach (var waternetline in Waternet.WaternetLineList)
                    {
                        if (waternetline.Name != "Surface line")
                        {
                            IList<GeometryPoint> intersectPoints =
                                surface.GeometrySurface.OuterLoop.ClosedGeometryIntersectionXzPointsWithGeometryPointList(waternetline.Points);
                            for (int i = 0; i < intersectPoints.Count; i++)
                            {
                                AddXcoordinate(intersectPoints[i].X);
                                changingPoints.Add(intersectPoints[i]);
                            }
                        }
                    }
                }
            }
        }

        /// <summary>
        /// This procedure determine the POP values at the given points (in table of Preconsolidation stresses) and the corresponding soil name
        /// and add it in popsInSoilArray.
        /// </summary>
        /// <param name="preconStress">Preconsolidation stress at a given point</param>
        private void DeterminePopFromPreconsolidationStressInSoil(PreConsolidationStress preconStress)
        {
            InterpolationColumn column = GetColumnAtX(preconStress.X);
            InterpolationArea area = column.getArea(preconStress.Z);
            if (area != null)
            {
                area.DetermineParameters();

                area.EffectiveStress = area.TotalStress - (area.PnLinePore + area.AanpasPore + area.PorePressureDueToDegreeofConsolidationLoad);

                var popCalculate = new CalculatePopSoil();
                popCalculate.Soil = area.Soil;
                popCalculate.XCoor = preconStress.X;
                popCalculate.ZCoor = preconStress.Z;
                popCalculate.Pop = area.DeterminePopFromPreconsolidationStress(preconStress.StressValue);

                popsInSoilArray.Add(popCalculate);
            }
            else
            {
                // The area can be null if the pre-consoldiation stress point is situated exactely on the surface line. 
                // To check this, the Z coordinate of the point is decreased by 0.0001 m, and a new search is performed.
                // If this new search still gives a null area, a log message is displayed.
                area = column.getArea(preconStress.Z - 0.0001);
                if (area != null)
                {
                    area.DetermineParameters();
                    area.EffectiveStress = area.TotalStress - (area.PnLinePore + area.AanpasPore + area.PorePressureDueToDegreeofConsolidationLoad);
                    var popCalculate = new CalculatePopSoil();
                    popCalculate.Soil = area.Soil;
                    popCalculate.XCoor = preconStress.X;
                    popCalculate.ZCoor = preconStress.Z;
                    popCalculate.Pop = area.DeterminePopFromPreconsolidationStress(preconStress.StressValue);
                    popsInSoilArray.Add(popCalculate);
                }
                else
                {
                    var logMessage = new LogMessage(LogMessageType.Warning, this, string.Format(LocalizationManager.GetTranslatedText(typeof(PreProcesStabilityCalculation),
                    "NotAbleToFindInWhichAreaTheFollowingPreconsolidationStressPointIsSituated")) + ": X = " + preconStress.X + ", Z = " + preconStress.Z);
                    if (LogManager.Messages.All(m => m.Message != logMessage.Message)) // Only add log message once
                    {
                        LogManager.Add(logMessage);
                    }
                }
            }
        }

        /// <summary>
        /// If a column has no area with a POP value (meaning no soil in this column where POP was calculated),
        /// the closest POP is used for those columns.
        /// </summary>
        private void FillNotUsedColumnsWithPop()
        {
            for (int k = 0; k < interpolationColumns.Count; k++)
            {
                if (interpolationColumns[k].Areas.Count > 0 && interpolationColumns[k].Areas[0].POP < -0.001)
                {
                    double xmid = 0.5*(interpolationColumns[k].XLeft + interpolationColumns[k].XRight);
                    if (Math.Abs(xmid - 41.77) < 0.1)
                    {
                        Console.WriteLine("slice");
                    }

                    double minDistance = 1000000.0;
                    double distance;
                    double pop = 0;
                    double number = 0;
                    for (int i = 0; i < popsInSoilArray.Count; i++)
                    {
                        distance = Math.Abs(xmid - popsInSoilArray[i].XCoor);
                        if (distance <= minDistance)
                        {
                            pop = pop + popsInSoilArray[i].Pop;
                            minDistance = distance;
                            number++;
                        }
                    }
                    pop = pop / number;
                    foreach (var area in interpolationColumns[k].Areas)
                    {
                        area.POP = pop;
                    }
                }
            }
        }

        /// <summary>
        /// Areas without POP (meaning no soil in column where POP was calculated) are filled by interpolation: 
        /// first by interpolation along the column if other areas with a POP value are present.
        /// If not, by horizontal interpolation, using the closest column. 
        /// </summary>
        private void DeterminePopForNotUsedSoilnames()
        {
            for (int k = 0; k < interpolationColumns.Count; k++)
            {
                FillNotUsedAreasInColumnsWithPopByInterpolation(interpolationColumns[k]);
            }

            FillNotUsedColumnsWithPop();
        }

        private void SetPopInAreasAtminusValue()
        {
            for (int i = 0; i < interpolationColumns.Count; i++)
            {
                foreach (var area in interpolationColumns[i].Areas)
                {
                    area.POP = -1;
                }
            }
        }

        /// <summary>
        /// First, this procedure determine the POP values at the given points (in table of Preconsolidation stresses). 
        /// Secondly, the areas having a soil name with is present in the table of Preconsolidation stresses are filled with the corresponding POP (in case
        /// a soil material has several POP values within a column, the average is used). Then, the other areas are filled with interpolation: by
        /// interpolation along the column if other areas with a POP value are present or by horizontal interpolation using the closest column if not.
        /// </summary>
        private void FillAreasWithPreconsolidationStressData()
        {
            // Intialize POP in all areas at -1
            SetPopInAreasAtminusValue();

            // Replace POP with POP from calculation
            for (int i = 0; i < preconStressArray.Count; i++)
            {
                // First fill all the POP in soil array with the corresponding soil name.
                // Note that a soil name can appear several times in popsInSoilArray.
                DeterminePopFromPreconsolidationStressInSoil(preconStressArray[i]);
            }

            var namesUsed = new List<Soil>(); // all the names where the pop is calculated
            for (int j = 0; j < popsInSoilArray.Count; j++)
            {
                Soil soil = popsInSoilArray[j].Soil;
                if (!namesUsed.Contains(soil))
                {
                    DeterminePopForNewSoilName(soil);
                    namesUsed.Add(soil);
                }
            }

            // now add pop to layers where no preconsolidation stress was calculated
            // this is done by interpolation, using the closest columns
            DeterminePopForNotUsedSoilnames();

            // for debug purpose ToDo Delete before final version 
//            for (int i = 0; i < interpolationColumns.Count; i++)
//            {
//                for (int j = 0; j < interpolationColumns[i].Areas.Count; j++)
//                {
//                    interpolationColumns[i].Areas[j].WriteData(writer);
//                }
//            }
        }

        /// <summary>
        /// If the checkbox "Use POP" is checked, the POP of the soil material given in the Materials table is used
        /// </summary>
        private void FillAreasInColumnsWithUserDefinedPop()
        {
            for (int k = 0; k < interpolationColumns.Count; k++)
            {
                for (int i = 0; i < interpolationColumns[k].Areas.Count; i++)
                {
                    if (interpolationColumns[k].Areas[i].Soil.UsePop)
                    {
                        interpolationColumns[k].Areas[i].POP = interpolationColumns[k].Areas[i].Soil.POP;
                    }
                }
            }
        }

        internal class CalculatePopSoil
        {
            public double Pop;
            public Soil Soil;
            public double XCoor;
            public double ZCoor;
        }
    }
}