//-----------------------------------------------------------------------
// <copyright file="StabilityCalculator.cs" company="Deltares">
//  Copyright (c) 2009 Deltares. All rights reserved.
// </copyright>
// <author>B.S.T.I.M. The</author>
// <email>tom.the@deltares.nl</email>
// <date>04-06-2010</date>
// <summary>Stability Calculator object</summary>
//-----------------------------------------------------------------------

using Deltares.Geotechnics;
using Deltares.MStab;
using Deltares.MStab.IO.Classic;
using Deltares.Stability;
using Deltares.Stability.Calculation;
using Deltares.Standard;
using Deltares.Standard.EventPublisher;
using Deltares.Standard.Language;
using System;
using System.Collections.Generic;
using System.Xml.Linq;
using System.IO;
using System.Reflection;
using System.Text.RegularExpressions;
using Deltares.Dam.Data.Assemblers;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Mathematics;
using Deltares.Soilbase;
using Deltares.Standard.IO.Xml;
using Deltares.Standard.Logging;

namespace Deltares.Dam.Data
{

    public class StabilityCalculator : IDisposable
    {
        protected ModelParametersForPLLines modelParametersForPLLines;
        private ProgramType programType;
        private double trafficLoad;
        private double minimalCircleDepth;
        private double requiredSafetyFactor;
        private string mstabProgramPath;
        private string slopeWProgramPath;
        private FailureMechanismeParamatersMStab failureMechanismeParamatersMStab;
        private IList<GaugePLLine> gaugePLLines;
        private IList<Gauge> gauges;
        private string stabilityProjectFilename = "";
        private List<LogMessage> errorMessages = new List<LogMessage>();

        public const string StabilitySubDir = @"Stability\";

        public StabilityCalculator(FailureMechanismeParamatersMStab failureMechanismeParamatersMStab,
            ProgramType programType,
            ModelParametersForPLLines modelParametersForPLLines, double trafficLoad, double minimalCircleDepth,
            double requiredSafetyFactor, string mstabProgramPath, string slopeWProgramPath,
            IList<GaugePLLine> gaugePLLines,
            IList<Gauge> gauges, SoilbaseDB soilBaseDB, SoilList soilList, ProbabilisticType probabilisticType)
        {
            this.programType = programType;
            this.modelParametersForPLLines = modelParametersForPLLines;
            this.trafficLoad = trafficLoad;
            this.minimalCircleDepth = minimalCircleDepth;
            this.requiredSafetyFactor = requiredSafetyFactor;
            this.mstabProgramPath = mstabProgramPath;
            this.slopeWProgramPath = slopeWProgramPath;
            // need to have own instance of failureMechanismeParamatersMStab because the calculation
            // makes changes to this object (e.g. failureMechanismeParamatersMStab.PLLines)
            // This will not function correctly when the calculations are performed
            // multi threaded (MWDAM-568)
            this.failureMechanismeParamatersMStab = failureMechanismeParamatersMStab.Clone();
            this.gaugePLLines = gaugePLLines;
            this.gauges = gauges;
            this.SoilBaseDB = soilBaseDB;
            this.SoilList = soilList;
            ProbabilisticType = probabilisticType;
            CalculationBaseDirectory = "";
        }

        public string CalculationBaseDirectory { get; set; }
        public SoilbaseDB SoilBaseDB { get; set; }
        public SoilList SoilList { get; set; }
        public ProbabilisticType ProbabilisticType { get; set; }
        public static string ModelSubDirectory = "";
        private MStabProject dotNetMstabProjectResults;
        public PhreaticAdaptionType? NWOPhreaticAdaption { get; set; }

        /// <summary>
        /// Create PLLines with selected model
        /// </summary>
        /// <param name="location">The location.</param>
        /// <param name="surfaceLine">The surface line.</param>
        /// <param name="soilProfile">The soil profile.</param>
        /// <param name="soilGeometry2DName">Name of the soil geometry2 d.</param>
        /// <param name="waterlevel">The waterlevel.</param>
        /// <param name="waterlevelLow">The waterlevel low.</param>
        /// <returns></returns>
        private PLLines CreateAllPLLines(Location location, SurfaceLine2 surfaceLine, SoilProfile soilProfile,
            string soilGeometry2DName, double waterlevel, double? waterlevelLow)
        {
            switch (location.PLLineCreationMethod)
            {
                case PLLineCreationMethod.ExpertKnowledgeLinearInDike:
                case PLLineCreationMethod.ExpertKnowledgeRRD:
                case PLLineCreationMethod.GaugesWithFallbackToExpertKnowledgeRRD:
                    return CreateAllPLLinesExpertKnowledge(location, surfaceLine, soilProfile, soilGeometry2DName, waterlevel, waterlevelLow);
                case PLLineCreationMethod.DupuitStatic:
                    return CreateAllPLLinesDupuit();
                case PLLineCreationMethod.DupuitDynamic:
                    throw new StabilityCalculationException("PL-Line creation with DupuitDynamic not yet implemented");
                default:
                    return null;
            }
        }

        /// <summary>
        /// Create PLLines with Dupuit model
        /// </summary>
        /// <returns></returns>
        private PLLines CreateAllPLLinesDupuit()
        {
            return null;
        }

        /// <summary>
        /// Create PLLines with expert knowledge
        /// </summary>
        /// <param name="location">The location.</param>
        /// <param name="surfaceLine">The surface line.</param>
        /// <param name="soilProfile">The soil profile.</param>
        /// <param name="soilGeometry2DName">Name of the soil geometry2 d.</param>
        /// <param name="waterlevel">The waterlevel.</param>
        /// <param name="waterlevelLow">The waterlevel low.</param>
        /// <returns></returns>
        private PLLines CreateAllPLLinesExpertKnowledge(Location location, SurfaceLine2 surfaceLine, SoilProfile soilProfile,
            string soilGeometry2DName, double waterlevel, double? waterlevelLow)
        {
            PLLines plLines = null;

            if (surfaceLine != null && surfaceLine.HasDike())
            {
                PLLinesCreator plLineCreator = new PLLinesCreator();
                plLineCreator.SurfaceLine = surfaceLine;
                plLineCreator.WaterLevelRiverHigh = waterlevel;
                plLineCreator.WaterLevelRiverLow = waterlevelLow;
                plLineCreator.IsUseLowWaterLevel = (failureMechanismeParamatersMStab.MStabParameters.GridPosition == MStabGridPosition.Left);
                plLineCreator.WaterLevelPolder = location.PolderLevel;
                plLineCreator.HeadInPLLine2 = location.HeadPL2;
                plLineCreator.HeadInPLLine3 = location.HeadPl3;
                plLineCreator.HeadInPLLine4 = location.HeadPl4;
                plLineCreator.ModelParametersForPLLines = this.modelParametersForPLLines;
                SoilProfile1D profile1D = soilProfile as SoilProfile1D;
                if (profile1D != null)
                {
                    plLineCreator.SoilProfile = profile1D;
                    plLineCreator.SoilGeometryType = SoilGeometryType.SoilGeometry1D;
                }
                else
                {
                    plLineCreator.SoilGeometryType = SoilGeometryType.SoilGeometry2D;
                    string mapForSoilGeometries2D = location.MapForSoilGeometries2D;
                    soilGeometry2DName = Path.Combine(Path.Combine(DamProject.ProjectMap, mapForSoilGeometries2D), soilGeometry2DName);
                    plLineCreator.SoilGeometry2DName = soilGeometry2DName;
                    plLineCreator.SoilBaseDB = this.SoilBaseDB;
                    plLineCreator.SoilList = this.SoilList;
                    plLineCreator.DikeEmbankmentMaterial =
                        this.SoilList.GetSoilByName(location.DikeEmbankmentMaterial);
                }
                plLineCreator.GaugePLLines = this.gaugePLLines;
                plLineCreator.Gauges = this.gauges;
                plLineCreator.IsAdjustPL3AndPL4SoNoUpliftWillOccurEnabled = true;
                // for stability this must be set to true
                plLineCreator.PlLineOffsetBelowDikeTopAtRiver = location.PlLineOffsetBelowDikeTopAtRiver;
                plLineCreator.PlLineOffsetBelowDikeTopAtPolder = location.PlLineOffsetBelowDikeTopAtPolder;
                plLineCreator.PlLineOffsetBelowShoulderBaseInside = location.PlLineOffsetBelowShoulderBaseInside;
                plLineCreator.PlLineOffsetBelowDikeToeAtPolder = location.PlLineOffsetBelowDikeToeAtPolder;
                plLineCreator.PlLineOffsetBelowDikeCrestMiddle = location.PlLineOffsetBelowDikeCrestMiddle;
                plLineCreator.PlLineOffsetFactorBelowShoulderCrest = location.PlLineOffsetFactorBelowShoulderCrest;
                plLineCreator.UsePlLineOffsetBelowDikeCrestMiddle = location.UsePlLineOffsetBelowDikeCrestMiddle;
                plLineCreator.UsePlLineOffsetFactorBelowShoulderCrest = location.UsePlLineOffsetFactorBelowShoulderCrest;
                plLineCreator.NWOPhreaticAdaption = NWOPhreaticAdaption;
                plLineCreator.XSoilGeometry2DOrigin = location.XSoilGeometry2DOrigin;

                plLines = plLineCreator.CreateAllPLLines(location);
            }

            return plLines;
        }

        public string GetStabilityCalculationDirectory()
        {
            string stabilityDirectory = GetStabilityCalculationBaseDirectory();
            if (StabilityCalculator.ModelSubDirectory != "")
                stabilityDirectory = Path.Combine(stabilityDirectory, StabilityCalculator.ModelSubDirectory);
            if (!Directory.Exists(stabilityDirectory))
                Directory.CreateDirectory(stabilityDirectory);
            return stabilityDirectory;
        }

        /// <summary>
        /// Get the directory where to create the MStab project files
        /// </summary>
        /// <returns></returns>
        private string GetStabilityCalculationBaseDirectory()
        {
            if (CalculationBaseDirectory == "")
            {
                CalculationBaseDirectory = DamProject.ProjectWorkingPath;
            }
            return Path.Combine(CalculationBaseDirectory, StabilitySubDir);
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="soilProfile"></param>
        /// <param name="soilGeometry2DName"></param>
        public void ConsistencyCheck(Scenario scenario, SoilProfile1D soilProfile, string soilGeometry2DName)
        {
            if (soilProfile != null)
            {
                if (soilProfile.BottomAquiferLayer == null)
                {
                    throw new DamFailureMechanismeCalculatorException(
                        String.Format("Soilprofile '{0}' does not contain aquifer layers.", soilProfile.Name));
                }
            }
        }

        private void AddForbiddenZoneToMstabProject(MStabProject mStabProject, Scenario scenario, SurfaceLine2 surfaceLine)
        {
            if (failureMechanismeParamatersMStab.MStabParameters.CalculationOptions.ZonesType.Equals(
                MStabZonesType.ForbiddenZone))
            {
                CreateForbiddenZone(scenario, surfaceLine);
                AddMissingSlipPlaneConstraintInfo(mStabProject.Stability.SlipPlaneConstraints);
            }
        }

        private void CalculateType(string stabilityProjectFilename, Scenario scenario, SoilProfile1D soilProfile,
            string soilGeometry2DName,
            StabilityServiceAgent stabilityServiceAgent, double riverLevel, MStabDesignEmbankment mstabDesignEmbankment)
        {

            //
            // Here start of implementation of DesignPoint Probabilistic Calculation
            // Perform DesignPoint Probabilistic Calculation when this.ProbabilisticType == ProbabilisticType.ProbabilisticAdvanced
            //
            if (dotNetMstabProjectResults != null)
            {
                dotNetMstabProjectResults.Dispose();
                dotNetMstabProjectResults = null;
            }
            failureMechanismeParamatersMStab.MStabParameters.CalculationOptions = new MStabCalculationOptions
            {
                MinimalCircleDepth = this.minimalCircleDepth,
                ZonesType = scenario.Location.StabilityZoneType
            };
            SoilProfile profile;
            var surfaceLine = scenario.GetMostRecentSurfaceLine(soilProfile,
                Path.GetFileName(soilGeometry2DName));
            //create mstabproj using .NET, without using delphi dll and xml files
            if (SelectedStabilityKernelType == StabilityKernelType.AdvancedWti)
            {
                profile = GetSoilProfileByType(soilProfile, scenario.Location.MapForSoilGeometries2D, soilGeometry2DName);
                UpdateProfileSoilsToProjectSoils(profile);

                var parameterValues = GetStabilityProjectParameterValues();
                var stabilityProjectObjectCreator = new StabilityProjectObjectCreator();  
                MStabProject mStabProject = stabilityProjectObjectCreator.CreateMacroStabilityProject(scenario, profile, failureMechanismeParamatersMStab, riverLevel, parameterValues);
                AddForbiddenZoneToMstabProject(mStabProject, scenario, surfaceLine);
                RunStabilityCalculationWtiKernel(mStabProject);
                return;
            }
            //use mstab stability kernel, create mstabproj using .NET
            if (SelectedStabilityKernelType == StabilityKernelType.AdvancedDotNet)
            {
                profile = GetSoilProfileByType(soilProfile, scenario.Location.MapForSoilGeometries2D, soilGeometry2DName);
                UpdateProfileSoilsToProjectSoils(profile);

                var parameterValues = GetStabilityProjectParameterValues();
                var stabilityProjectObjectCreator = new StabilityProjectObjectCreator();
                MStabProject mStabProject = stabilityProjectObjectCreator.CreateMacroStabilityProject(scenario,
                    profile, failureMechanismeParamatersMStab, riverLevel, parameterValues);
                AddForbiddenZoneToMstabProject(mStabProject, scenario, surfaceLine);
                var projectString = SaveStabilityProjectFileToString(mStabProject);

                RunMstabStabilityCalculation(projectString);
                return;
            }

            XDocument mstabXML = CreateMStabXmlDoc(stabilityProjectFilename, scenario, soilProfile, soilGeometry2DName,
                riverLevel, mstabDesignEmbankment, surfaceLine);
            mstabXML.Save(stabilityProjectFilename + ".xml");

            stabilityServiceAgent.CreateProjectFile(mstabXML.ToString());
            // use the current stability kernel
            if (SelectedStabilityKernelType == StabilityKernelType.DamClassic)
            {
                if (programType == ProgramType.SlopeW)
                {
                    stabilityServiceAgent.CalculateSlopeWProject(stabilityProjectFilename);
                }
                else
                {
                    stabilityServiceAgent.CalculateMStabProject(stabilityProjectFilename);
                }
            }
            // use the converted .NET stability kernel 
            else if (SelectedStabilityKernelType == StabilityKernelType.DamClassicDotNet)
            {
                RunMstabStabilityCalculation(stabilityProjectFilename);
            }

            // use the WTI Stability Kernel, create mstabproj using delphi
            else if (SelectedStabilityKernelType == StabilityKernelType.DamClassicWti)
            {
                MStabProject mStabProject = ReadStabilityModel(stabilityProjectFilename);
                mStabProject.Stability.Location.WaterLevelRiver = riverLevel;
                AddMissingSlipPlaneConstraintInfo(mStabProject.Stability.SlipPlaneConstraints);
                mStabProject.Stability.SurfaceLine2 = surfaceLine;
                RunStabilityCalculationWtiKernel(mStabProject);
            }
        }

        /// <summary>
        /// Updates the profile soils to project soils.
        /// </summary>
        /// <param name="profile">The profile.</param>
        /// <exception cref="StabilityCalculationException"></exception>
        private void UpdateProfileSoilsToProjectSoils(SoilProfile profile)
        {
            // Assign the soils from the soillist to the soils in the profile
            var soilProfile2D = profile as SoilProfile2D;
            if (soilProfile2D != null)
            {
                foreach (var surface in soilProfile2D.Surfaces)
                {
                    var soil = SoilList.GetSoilByName(surface.Soil.Name);
                    if (soil == null)
                    {
                        throw new StabilityCalculationException(String.Format("No matching soil found for {0}", surface.Soil.Name));
                    }
                    surface.Soil = soil;
                }
            }
            else
            {
                var soilProfile1D = profile as SoilProfile1D;
                if (soilProfile1D != null)
                {
                    foreach (var layer in soilProfile1D.Layers)
                    {
                        var soil = SoilList.GetSoilByName(layer.Soil.Name);
                        if (soil == null)
                        {
                            throw new StabilityCalculationException(String.Format("No matching soil found for {0}", layer.Soil.Name));
                        }
                        layer.Soil = soil;
                    }
                }
                else
                {
                    throw new StabilityCalculationException(String.Format("Profile '{0}' has unexpected soilprofile type ", profile.Name));
                }
            }
        }

        private void AddMissingSlipPlaneConstraintInfo(SlipplaneConstraints slipPlaneConstraints)
        {
            slipPlaneConstraints.CreateZones =
                failureMechanismeParamatersMStab.MStabParameters.ZonesType == MStabZonesType.ForbiddenZone;
            if (slipPlaneConstraints.CreateZones)
            {
                slipPlaneConstraints.SlipPlaneMinDepth = minimalCircleDepth;

                if (failureMechanismeParamatersMStab.MStabParameters.ForbiddenZone.IsXEntryMaxUsed)
                {
                    slipPlaneConstraints.XEntryMax = failureMechanismeParamatersMStab.MStabParameters.ForbiddenZone.XEntryMax;
                }
                else
                {
                    slipPlaneConstraints.XEntryMax = double.NaN;
                }

                if (failureMechanismeParamatersMStab.MStabParameters.ForbiddenZone.IsXEntryMinUsed)
                {
                    slipPlaneConstraints.XEntryMin = failureMechanismeParamatersMStab.MStabParameters.ForbiddenZone.XEntryMin;
                }
                else
                {
                    slipPlaneConstraints.XEntryMin = double.NaN;
                }
            }
        }

        private static SoilProfile GetSoilProfileByType(SoilProfile1D soilProfile, string mapForSoilGeometries2D, string soilGeometry2DName)
        {
            SoilProfile profile = null;
            if (soilProfile != null)
            {
                profile = soilProfile;
            }
            else
            {
                var projDir = Path.GetDirectoryName(DamProject.ProjectWorkingPath);
                var fullPath = Path.Combine(projDir, mapForSoilGeometries2D, Path.GetFileName(soilGeometry2DName));
                // Bka: for now, see if there is a dsx version too and if so prefer that.
                var dsxPath = Path.ChangeExtension(fullPath, ".dsx");
                MStabProject mStab = null;
                //if (soilGeometry2DName.EndsWith(".dsx", StringComparison.OrdinalIgnoreCase))
                if (File.Exists(dsxPath))
                {
                    try
                    {
                        var xmlDeserializer = new XmlDeserializer();
                        mStab = (MStabProject) xmlDeserializer.XmlDeserialize(dsxPath, typeof (MStabProject));
                    }
                    catch (Exception ex)
                    {
                        LogManager.Add(new LogMessage(LogMessageType.Error, typeof (Converter), ex.Message));
                    }
                    if (mStab != null)
                    {
                        profile = mStab.Stability.SoilProfile;
                        mStab.Dispose();
                    }
                }
                else
                {
                    var converter = new Converter();
                    mStab = new MStabProject();
                    mStab.Stability = converter.ConvertClassicMStab(fullPath);
                    if (mStab.Stability.SoilProfile.Surfaces != null && mStab.Stability.SoilProfile.Surfaces.Count > 0)
                    {
                        // Ensure that the bottom layer of the geometry is seen as Aquifer.
                        // There is no other way to get aquifer information as Delphi DGS does not work with aquifers and also not 
                        // with soiltypes so you can not derive anything from that too.
                        mStab.Stability.SoilProfile.Surfaces[0].IsAquifer = true;
                    }
                    profile = mStab.Stability.SoilProfile;
                    mStab.Dispose();
                }
            }
            return profile;
        }

        /// <summary>
        /// get stability project input parameters from DAM input
        /// </summary>
        /// <returns></returns>
        private StabilityProjectParameterValues GetStabilityProjectParameterValues()
        {
            StabilityProjectParameterValues parameterValues = new StabilityProjectParameterValues();
            parameterValues.MinimalCircleDepth = minimalCircleDepth;
            parameterValues.TrafficLoad = trafficLoad;
            parameterValues.PhreaticAdaptionType = NWOPhreaticAdaption != null
                ? (PhreaticAdaptionType)NWOPhreaticAdaption
                : PhreaticAdaptionType.None;
            parameterValues.PLlineCreationMethod = modelParametersForPLLines.PLLineCreationMethod;
            return parameterValues;
        }

        /// <summary>
        /// run stability calculation using WTI stability kernel
        /// </summary>
        /// <param name="mStabProject"></param>
        private void RunStabilityCalculationWtiKernel(MStabProject mStabProject)
        {
            StabilityModel stabilityModel = mStabProject.Stability;
            stabilityModel.CalculationModel = CalculationModel.RTO;
            var calculation = new Stability.Calculation2.StabilityCalculation();
            calculation.RunCalculation(stabilityModel);
            // Todo mstabproject still to be filled with the results from the stabilityModel calculation.
            if (stabilityModel.MinimumSafetyCurve != null)
            {
                mStabProject.Result = CalculationResult.Succeeded;
                mStabProject.SlidingData = new SlidingModel() { CurrentZone = new Zone() { MinimumSafetyCurve = stabilityModel.MinimumSafetyCurve } };

            }
            dotNetMstabProjectResults = mStabProject;
        }

        /// <summary>
        /// run stability calculation using classic .NET kernel 
        /// </summary>
        /// <param name="projectString"></param>
        private void RunMstabStabilityCalculation(string projectString)
        {
            using (var calculation = new Stability.Calculation.StabilityCalculation())
            {
                CalculationResult loadResult = calculation.Load(projectString);
                MStabProject loadedMStabProject = calculation.ClaimLoadedMStabProject();
                if (loadResult == CalculationResult.Succeeded)
                {
                    StabilityModel stabilityModel = loadedMStabProject.Stability;
                    stabilityModel.CalculationModel = CalculationModel.DSerie;                    
                    CalculationResult runResult = calculation.Run();
                    if (runResult == CalculationResult.Succeeded)
                    {
                        string results = null;
                        CalculationResult calculationResult = calculation.GetResults(ref results);
                        if (calculationResult == CalculationResult.Succeeded)
                        {
                            loadedMStabProject.ProcessResults(results);
                            dotNetMstabProjectResults = loadedMStabProject;
                        }
                        else
                        {
                            loadedMStabProject.Dispose();
                            throw new DamFailureMechanismeCalculatorException(calculation.Messages[0].Message);
                        }
                    }
                }
            }
        }

        public MStabProject GetDotNetStabilityKernelResults()
        {
            return dotNetMstabProjectResults;
        }

        private MStabProject ReadStabilityModel(string fileName)
        {
            MStabProject project = null;
            DataEventPublisher.InvokeWithoutPublishingEvents(() =>
            {
                if (!fileName.EndsWith(".dsx"))
                {
                    var converter = new Converter();
                    project = new MStabProject();
                    project.Stability = converter.ConvertClassicMStab(fileName);
                }
                else
                {
                    XmlDeserializer deserializer = new XmlDeserializer();
                    project = (MStabProject)deserializer.XmlDeserialize(fileName, typeof(MStabProject));
                }
            });
            
            return project;
        }

        /// <summary>
        /// Determine a filename based on the different inout parameters
        /// </summary>
        /// <param name="locationName"></param>
        /// <param name="scenarioName"></param>
        /// <param name="soilGeometryName"></param>
        /// <param name="iterationIndex"></param>
        /// <returns></returns>
        public static string DetermineCalculationFilename(string locationName, string scenarioName,
            string soilGeometryName, int iterationIndex)
        {
            string calculationName;
            if (iterationIndex < 0)
            {
                calculationName = String.Format("Loc({0})_Sce({1})_Pro({2})", locationName, scenarioName,
                    soilGeometryName);
            }
            else
            {
                calculationName = String.Format("Loc({0})_Sce({1})_Pro({2})_Ite({3})", locationName, scenarioName,
                    soilGeometryName, iterationIndex);
            }
            return Regex.Replace(calculationName, @"[\\\/:\*\?""'<>|.]", "_");
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="soilProfile"></param>
        /// <param name="soilGeometry2DName"></param>
        /// <param name="iterationIndex"></param>
        public void Calculate(Scenario scenario, SoilProfile1D soilProfile, string soilGeometry2DName,
            int iterationIndex)
        {
            Calculate(scenario, soilProfile, soilGeometry2DName, iterationIndex, null);
        }

        /// <summary>
        /// Calculate 1 scenario 
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="soilProfile"></param>
        /// <param name="soilGeometry2DName"></param>
        /// <param name="iterationIndex"></param>
        public void Calculate(Scenario scenario, SoilProfile1D soilProfile, string soilGeometry2DName,
            int iterationIndex, MStabDesignEmbankment mstabDesignEmbankment)
        {
            failureMechanismeParamatersMStab.MStabParameters.IsProbabilistic = ((this.ProbabilisticType ==
                                                                                 ProbabilisticType.Probabilistic) ||
                                                                                (this.ProbabilisticType ==
                                                                                 ProbabilisticType
                                                                                     .ProbabilisticFragility));
            string soilGeometryName;
            if (!string.IsNullOrEmpty(soilGeometry2DName))
            {
                soilGeometryName = Path.GetFileName(soilGeometry2DName);
            }
            else
            {
                soilGeometryName = soilProfile.Name;
            }

            string calculationName = DetermineCalculationFilename(scenario.Location.Name, scenario.LocationScenarioID,
                soilGeometryName, iterationIndex);

            try
            {
                StabilityServiceAgent stabilityServiceAgent = new StabilityServiceAgent();
                stabilityServiceAgent.ProgramType = this.programType;
                stabilityServiceAgent.SlopeWExePath = this.slopeWProgramPath;
                stabilityServiceAgent.MStabExePath = this.mstabProgramPath;
                string stabilityDirectory = GetStabilityCalculationDirectory();
                string filenameExtension = GetFilenameExtension();
                string fileName = calculationName + filenameExtension;
                stabilityProjectFilename = Path.Combine(stabilityDirectory, fileName);

                MStabResults mStabResults = new MStabResults();
                mStabResults.Init();
                switch (ProbabilisticType)
                {
                    case ProbabilisticType.Deterministic:
                        {
                            CalculateType(stabilityProjectFilename, scenario, soilProfile, soilGeometry2DName,
                                stabilityServiceAgent, scenario.RiverLevel, mstabDesignEmbankment);
                            scenario.SetFailureProbabilityStability(soilProfile, Path.GetFileName(soilGeometry2DName), null);
                            break;
                        }
                    case ProbabilisticType.Probabilistic:
                        {
                            CalculateType(stabilityProjectFilename, scenario, soilProfile, soilGeometry2DName,
                                stabilityServiceAgent, scenario.RiverLevel, mstabDesignEmbankment);
                            double beta = 0;
                            beta = stabilityServiceAgent.ExtractBeta(stabilityProjectFilename);
                            scenario.SetFailureProbabilityStability(soilProfile, soilGeometry2DName,
                                Probabilistic.Probabilistic.NormalDistribution(-beta));
                            break;
                        }
                    case ProbabilisticType.ProbabilisticFragility:
                        {
                            double[] waterLevels = scenario.DetermineProperWaterlevelsForProbabilisticAdvanced();
                            double[] betas = new double[3];
                            for (int i = 0; i < 3; i++)
                            {
                                string fileNameAdvanced = calculationName + "_WL" + i + ".sti";
                                stabilityProjectFilename = Path.Combine(stabilityDirectory, fileNameAdvanced);
                                CalculateType(stabilityProjectFilename, scenario, soilProfile, soilGeometry2DName,
                                    stabilityServiceAgent, waterLevels[i], mstabDesignEmbankment);
                                betas[i] = stabilityServiceAgent.ExtractBeta(stabilityProjectFilename);
                            }
                            var designPointWater = new DesignPointCalculation();
                            designPointWater.Waterlevels = waterLevels;
                            designPointWater.Betas = betas;
                            // Note Bka: For now, set MHW to original max water level (= river level + WaterHeightDecimeringsHoogte) and set
                            //           Decimate to the original WaterHeightDecimeringsHoogte. Han Best has to approve this! 
                            designPointWater.MHW = (double)(scenario.RiverLevel + scenario.WaterHeightDecimeringsHoogte);
                            designPointWater.Decimate = (double)scenario.WaterHeightDecimeringsHoogte;
                            designPointWater.Exceed = DesignPointCalculation.ExceedingSet.twoThousend;
                            designPointWater.IsMaxLevelUsed = false;
                            designPointWater.MaxLevel = 0;

                            if (designPointWater.CalculateTheWaterDesignpoint())
                            {
                                scenario.SetFailureProbabilityStability(soilProfile, soilGeometry2DName,
                                    Probabilistic.Probabilistic.NormalDistribution(-designPointWater.Beta));
                            }
                            break;
                        }
                }
                if (SelectedStabilityKernelType == StabilityKernelType.DamClassic)
                {
                    mStabResults = stabilityServiceAgent.ExtractStabilityResults(stabilityProjectFilename);
                    mStabResults = SetMStabAdministrationResults(mStabResults, iterationIndex, calculationName);
                    DetermineMStabResultsEntryPoint(ref mStabResults,
                        failureMechanismeParamatersMStab.MStabParameters.GridPosition);
                    DetermineMStabResultsExitPoint(ref mStabResults,
                        failureMechanismeParamatersMStab.MStabParameters.GridPosition);
                    scenario.SetMStabResults(soilProfile, Path.GetFileName(soilGeometry2DName), mStabResults);
                }
                else
                {
                    if (dotNetMstabProjectResults != null)
                    {
                        string stabilityResultsFileName = stabilityProjectFilename.Replace("sti", "dsx");
                        if (File.Exists(stabilityResultsFileName))
                        {
                            File.Delete(stabilityResultsFileName);
                        }
                        var filledGeometry = dotNetMstabProjectResults.Geometry;
                        // Todo het juist vullen van het Mstab project met alle waarden op de benodige plek (iom Rob) dus ook de results                        
                        SaveStabilityProjectToDsxFile(stabilityResultsFileName, dotNetMstabProjectResults);
                        mStabResults = ExtractDotNetStabilityKernelResults();
                        mStabResults = SetMStabAdministrationResults(mStabResults, iterationIndex, calculationName);
                        if (dotNetMstabProjectResults.Stability.MinimumSafetyCurve == null &&
                            !IsSlidingDataMinimumSafetyCurveAvailable())
                        {
                            return;
                        }
                        var cn = Path.GetFileNameWithoutExtension(stabilityResultsFileName);
                        mStabResults.CalculationName = cn;
                        DetermineDotNetStabilityKernelResultsEntryPoint(mStabResults);
                        // TODO: Exit point is not read correctly
                        DetermineDotNetStabilityKernelResultsExitPoint(mStabResults);
                        scenario.SetMStabResults(soilProfile, Path.GetFileName(soilGeometry2DName), mStabResults);
                    }
                }
            }

            catch (Exception e)
            {
                throw new DamFailureMechanismeCalculatorException(
                    String.Format("Error calculating stability factor for '{0}'", calculationName), e);
            }
        }

        private static MStabResults SetMStabAdministrationResults(MStabResults mStabResults, int iterationIndex,
            string calculationName)
        {
            mStabResults.CalculationName = calculationName;
            mStabResults.CalculationSubDir = StabilitySubDir;
            mStabResults.IterationNumber = iterationIndex;
            if (StabilityCalculator.ModelSubDirectory != "")
            {
                mStabResults.CalculationSubDir = Path.Combine(mStabResults.CalculationSubDir,
                    StabilityCalculator.ModelSubDirectory);
            }
            return mStabResults;
        }


        private bool IsSlidingDataMinimumSafetyCurveAvailable()
        {
            var notAvailable = dotNetMstabProjectResults.SlidingData == null ||
                   dotNetMstabProjectResults.SlidingData.CurrentZone == null ||
                   dotNetMstabProjectResults.SlidingData.CurrentZone.MinimumSafetyCurve == null;
            return !notAvailable;
        }

        private MStabResults ExtractDotNetStabilityKernelResults()
        {
            MStabResults mStabResults = new MStabResults();
            if (dotNetMstabProjectResults.Stability.HasZonePlot)
            {
                throw new NotImplementedException("Zoning features are not implementen yet in the Wti Stability kernel");
            }
            else
            {
                if (dotNetMstabProjectResults.Stability.MinimumSafetyCurve == null)
                {
                    if (!IsSlidingDataMinimumSafetyCurveAvailable())
                    {
                        return mStabResults;
                    }
                    else
                    {
                        dotNetMstabProjectResults.Stability.MinimumSafetyCurve =
                            dotNetMstabProjectResults.SlidingData.CurrentZone.MinimumSafetyCurve;
                    }
                }
                double safetyFactor = dotNetMstabProjectResults.Stability.MinimumSafetyCurve.SafetyFactor;
                mStabResults.zone1.safetyFactor = safetyFactor;
                mStabResults.zone1.circleSurfacePointLeftXCoordinate = dotNetMstabProjectResults.Stability.MinimumSafetyCurve.LeftPoint.X;
                mStabResults.zone1.circleSurfacePointRightXCoordinate = dotNetMstabProjectResults.Stability.MinimumSafetyCurve.RightPoint.X;
                if (dotNetMstabProjectResults.Stability.ModelOption == ModelOptions.UpliftVan ||
                    dotNetMstabProjectResults.Stability.ModelOption == ModelOptions.UpliftSpencer)
                {
                    mStabResults.zone1.safetyFactor = safetyFactor / 1.05;
                }
            }
            return mStabResults;
        }

        private void DetermineDotNetStabilityKernelResultsExitPoint(MStabResults results)
        {
            if (dotNetMstabProjectResults.Stability.GridOrientation == GridOrientation.Outwards)
            {
                if (dotNetMstabProjectResults.Stability.MinimumSafetyCurve.LeftPoint != null)
                {
                    // Outward stability
                    results.zone1.exitPointXCoordinate = dotNetMstabProjectResults.Stability.MinimumSafetyCurve.LeftPoint.X;
                    if (results.zone2.HasValue)
                    {
                    }
                }
            }
            else
            {
                // Inward stability
                if (dotNetMstabProjectResults.Stability.MinimumSafetyCurve.RightPoint != null)
                {
                    results.zone1.exitPointXCoordinate = dotNetMstabProjectResults.Stability.MinimumSafetyCurve.RightPoint.X;
                    if (results.zone2.HasValue)
                    {
                    }
                }
            }
        }

        private void DetermineDotNetStabilityKernelResultsEntryPoint(MStabResults mStabResults)
        {
            if (dotNetMstabProjectResults.Stability.GridOrientation == GridOrientation.Outwards)
            {
                if (dotNetMstabProjectResults.Stability.MinimumSafetyCurve.RightPoint != null)
                {
                    mStabResults.zone1.entryPointXCoordinate =
                        dotNetMstabProjectResults.Stability.MinimumSafetyCurve.RightPoint.X;
                    if (mStabResults.zone2.HasValue)
                    {
                        /*MStabResultsSingleZone zone = mstabResults.zone2.Value;
                            zone.entryPointXCoordinate = zone.circleSurfacePointRightXCoordinate;
                            mstabResults.zone2 = zone;*/
                    }
                }
            }
            else
            {
                if (dotNetMstabProjectResults.Stability.MinimumSafetyCurve.LeftPoint != null)
                {
                    // Inward stability
                    mStabResults.zone1.entryPointXCoordinate =
                        dotNetMstabProjectResults.Stability.MinimumSafetyCurve.LeftPoint.X;
                    if (mStabResults.zone2.HasValue)
                    {
                        /*MStabResultsSingleZone zone = mstabResults.zone2.Value;
                            zone.entryPointXCoordinate = zone.circleSurfacePointLeftXCoordinate;
                            mstabResults.zone2 = zone;*/
                    }
                }
            }
        }

        /// <summary>
        /// Determine filename extension according to which failuremechanism is used
        /// </summary>
        /// <returns></returns>
        public string GetFilenameExtension()
        {
            string filenameExtension;
            if (programType == ProgramType.SlopeW)
            {
                filenameExtension = ".xml";
            }
            else
            {
                filenameExtension = ".sti";
            }
            return filenameExtension;
        }

        /// <summary>
        /// Depending on outward/inward stability (position of grid) the entrypoint of the slipcricle is determined
        /// </summary>
        /// <param name="mstabResults"></param>
        /// <param name="mstabGridPosition"></param>
        private void DetermineMStabResultsEntryPoint(ref MStabResults mstabResults, MStabGridPosition mstabGridPosition)
        {
            if (mstabGridPosition == MStabGridPosition.Left)
            {
                // Outward stability
                mstabResults.zone1.entryPointXCoordinate = mstabResults.zone1.circleSurfacePointRightXCoordinate;
                if (mstabResults.zone2.HasValue)
                {
                    MStabResultsSingleZone zone = mstabResults.zone2.Value;
                    zone.entryPointXCoordinate = zone.circleSurfacePointRightXCoordinate;
                    mstabResults.zone2 = zone;
                }
            }
            else
            {
                // Inward stability
                mstabResults.zone1.entryPointXCoordinate = mstabResults.zone1.circleSurfacePointLeftXCoordinate;
                if (mstabResults.zone2.HasValue)
                {
                    MStabResultsSingleZone zone = mstabResults.zone2.Value;
                    zone.entryPointXCoordinate = zone.circleSurfacePointLeftXCoordinate;
                    mstabResults.zone2 = zone;
                }
            }
        }


        /// <summary>
        /// Depending on outward/inward stability (position of grid) the exitpoint of the slipcricle is determined
        /// </summary>
        /// <param name="mstabResults"></param>
        /// <param name="mstabGridPosition"></param>
        private void DetermineMStabResultsExitPoint(ref MStabResults mstabResults, MStabGridPosition mstabGridPosition)
        {
            if (mstabGridPosition == MStabGridPosition.Left)
            {
                // Outward stability
                mstabResults.zone1.exitPointXCoordinate = mstabResults.zone1.circleSurfacePointLeftXCoordinate;
                if (mstabResults.zone2.HasValue)
                {
                    MStabResultsSingleZone zone = mstabResults.zone2.Value;
                    zone.exitPointXCoordinate = zone.circleSurfacePointLeftXCoordinate;
                    mstabResults.zone2 = zone;
                }
            }
            else
            {
                // Inward stability
                mstabResults.zone1.exitPointXCoordinate = mstabResults.zone1.circleSurfacePointRightXCoordinate;
                if (mstabResults.zone2.HasValue)
                {
                    MStabResultsSingleZone zone = mstabResults.zone2.Value;
                    zone.exitPointXCoordinate = zone.circleSurfacePointRightXCoordinate;
                    mstabResults.zone2 = zone;
                }
            }
        }

        private void CreateForbiddenZone(Scenario scenario, SurfaceLine2 surfaceLine)
        {
            var dikeTopAtPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            // Zonestype is ForbiddenZone; TODO: Combine with code in StabilityCalculator?
            double maxZoneX = dikeTopAtPolder.X +
                              scenario.Location.ForbiddenZoneFactor * (surfaceLine.GetDikeToeInward().X - dikeTopAtPolder.X);
            failureMechanismeParamatersMStab.MStabParameters.ForbiddenZone = new MStabForbiddenZone
            {
                IsXEntryMinUsed = false,
                XEntryMin = 0.0,
                IsXEntryMaxUsed = true,
                XEntryMax = maxZoneX
            };                     
        }

        /// <summary>
        /// Create XML definition for Stability calculation
        /// </summary>
        /// <param name="mstabProjectFilename"></param>
        /// <param name="scenario"></param>
        /// <param name="soilProfile"></param>
        /// <param name="soilGeometry2DName"></param>
        /// <param name="riverLevel"></param>
        /// <param name="mstabDesignEmbankment"></param>
        /// <param name="surfaceLine"></param>
        /// <returns></returns>
        public XDocument CreateMStabXmlDoc(string mstabProjectFilename, Scenario scenario, SoilProfile soilProfile,
            string soilGeometry2DName, double riverLevel,
            MStabDesignEmbankment mstabDesignEmbankment, SurfaceLine2 surfaceLine)
        {
            SoilProfile1D profile1D = soilProfile as SoilProfile1D;
            SoilProfile2D profile2D = soilProfile as SoilProfile2D;
            ConsistencyCheck(scenario, profile1D, soilGeometry2DName);
            failureMechanismeParamatersMStab.Location = scenario.Location;
            if (profile1D != null)
            {
                failureMechanismeParamatersMStab.SoilProfile = profile1D;
                // 1d-geometry
                failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.SoilGeometryType =
                    SoilGeometryType.SoilGeometry1D;
            }
            else
            {
                // 2d-geometry
                failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.SoilGeometryType =
                    SoilGeometryType.SoilGeometry2D;
            }
            // Geometry Creation Options
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.SoilGeometry2DFilename =
                soilGeometry2DName;

            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.MaterialForDike =
                scenario.Location.DikeEmbankmentMaterial;
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.MaterialForShoulder =
                scenario.Location.ShoulderEmbankmentMaterial;
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.IsUseOriginalPLLineAssignments =
                scenario.Location.IsUseOriginalPLLineAssignments;
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.IsDesign =
                (mstabDesignEmbankment != null);
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.XOffsetSoilGeometry2DOrigin =
                -scenario.Location.XSoilGeometry2DOrigin;
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.PLLineAssignment = 
                CalculationHelper.PLLineCreationMethod2PLLineAssignment(scenario.Location.PLLineCreationMethod);
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.IntrusionVerticalWaterPressureType =
                scenario.Location.IntrusionVerticalWaterPressure.Value;
            failureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.PenetrationLength =
                scenario.Location.PenetrationLength;
            // End of Geometry Creation Options
            // Design options
            failureMechanismeParamatersMStab.Design = mstabDesignEmbankment;

            failureMechanismeParamatersMStab.SurfaceLine = surfaceLine;
            failureMechanismeParamatersMStab.RiverLevel = riverLevel; // scenario.RiverLevel;
            failureMechanismeParamatersMStab.DikeTableHeight = 
                scenario.DikeTableHeight ?? surfaceLine.GetDefaultDikeTableHeight() ?? 0;
            failureMechanismeParamatersMStab.TrafficLoad = this.trafficLoad;

            // Horizontal balance; TODO: Combine with code in StabilityCalculation
            if (failureMechanismeParamatersMStab.MStabParameters.Model == MStabModelType.HorizontalBalance)
            {
                if (profile1D == null)
                {
                    throw new DamFailureMechanismeCalculatorException(
                        "Model horizontal balance does not support 2d-geometries");
                }
                failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea = new HorizontalBalanceArea();
                failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.XLeft =
                    surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X;
                failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.XRight =
                    surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
                failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.YTop = riverLevel;
                failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.YBottom =
                    profile1D.InBetweenAquiferLayer != null
                        ? profile1D.InBetweenAquiferLayer.TopLevel
                        : profile1D.BottomAquiferLayer.TopLevel;
                int planeCount =
                    (int)(Math.Round((failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.YTop -
                                      failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.YBottom) / 0.25));
                failureMechanismeParamatersMStab.MStabParameters.HorizontalBalanceArea.PlaneCount = Math.Min(planeCount, 50);

            }

            // Zonestype is ZoneAreas; TODO: Combine with code in StabilityCalculation
            var dikeTopAtPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
            if (
                failureMechanismeParamatersMStab.MStabParameters.CalculationOptions.ZonesType.Equals(
                    MStabZonesType.ZoneAreas))
            {
                double? dikeTableHeight = scenario.DikeTableHeight ?? surfaceLine.GetDefaultDikeTableHeight() ?? null;
                if (!dikeTableHeight.HasValue)
                    throw new DamFailureMechanismeCalculatorException("Surface line has no dike table height.");
                failureMechanismeParamatersMStab.MStabParameters.ZoneAreas = new MStabZoneAreas
                {
                    DikeTableHeight = dikeTableHeight.Value,
                    DikeTableWidth = scenario.Location.ZoneAreaRestSlopeCrestWidth,
                    SafetyFactorZone1A = scenario.ModelFactors.RequiredSafetyFactorStabilityInnerSlope ?? this.requiredSafetyFactor,
                    SafetyFactorZone1B = scenario.ModelFactors.RequiredSafetyFactorStabilityInnerSlope ?? this.requiredSafetyFactor,
                    XCoordinateDikeTopAtPolder = dikeTopAtPolder.X,
                    XCoordinateDikeTopAtRiver = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X,
                    XCoordinateStartRestProfile = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).X
                };
            }

            if (failureMechanismeParamatersMStab.MStabParameters.CalculationOptions.ZonesType.Equals(
                MStabZonesType.ForbiddenZone))
            {
                CreateForbiddenZone(scenario, surfaceLine);
            }

            // Make sure riverlevel is correct with respect to surfaceline
            double riverLevelLow = double.NaN;
            if (failureMechanismeParamatersMStab.MStabParameters.GridPosition == MStabGridPosition.Left && scenario.RiverLevelLow.HasValue)
            {
                riverLevelLow = scenario.RiverLevelLow.Value;
            }
            double? riverLevelHigh = riverLevel;
            var surfaceLevelOutside = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside);
            if (riverLevelHigh < surfaceLevelOutside.Z)
            {
                var riverLevelHighIsBelowSurfaceLevelOutside = Path.GetFileName(mstabProjectFilename) + ": " +
                                                               LocalizationManager.GetTranslatedText(this.GetType(),
                                                                   "riverLevelHighIsBelowSurfaceLevelOutside");
                LogMessage logMessage = new LogMessage(LogMessageType.Warning, null,
                    String.Format(riverLevelHighIsBelowSurfaceLevelOutside, riverLevelHigh,
                        surfaceLevelOutside.Z));
                errorMessages.Add(logMessage);
            }
            
            var currentSurfaceLine = scenario.GetMostRecentSurfaceLine(soilProfile, Path.GetFileName(soilGeometry2DName));
            if (currentSurfaceLine == null)
            {
                currentSurfaceLine = surfaceLine;
            }
            if (SelectedStabilityKernelType == StabilityKernelType.AdvancedWti || 
                SelectedStabilityKernelType == StabilityKernelType.AdvancedDotNet)
            {
                var lphreaticAdaptionType = NWOPhreaticAdaption != null
                ? (PhreaticAdaptionType)NWOPhreaticAdaption
                : PhreaticAdaptionType.None;
                var stabilityProjectObjectCreator = new StabilityProjectObjectCreator();
                failureMechanismeParamatersMStab.PLLines = stabilityProjectObjectCreator.CreateAllPlLinesUsingWaternetCreator(scenario.Location,
                    currentSurfaceLine, soilProfile, lphreaticAdaptionType, modelParametersForPLLines.PenetrationLength, riverLevelHigh.Value, modelParametersForPLLines.PLLineCreationMethod, riverLevelLow);
            }
            else
            {
                failureMechanismeParamatersMStab.PLLines = CreateAllPLLines(scenario.Location, currentSurfaceLine,
                    soilProfile, soilGeometry2DName, riverLevelHigh.Value, riverLevelLow);
            }
            // Slip circle definition for Uplift Van; TODO: Combine with code in StabilityCalculation
            if (this.failureMechanismeParamatersMStab.MStabParameters.Model == MStabModelType.UpliftVan)
            {
                // Determine right side of slip plane grid (right grid)
                // This is the location with the lowest uplift factor or, if present, the second NWO point
                var upliftLocationAndResult = this.GetLocationWithLowestUpliftFactor(currentSurfaceLine, soilProfile,
                    soilGeometry2DName,
                    failureMechanismeParamatersMStab.PLLines,
                    scenario.Location);
                double upliftCriterion =
                    scenario.GetUpliftCriterionStability(scenario.Location.ModelFactors.UpliftCriterionStability);
                bool isUplift = !(upliftLocationAndResult == null) &&
                                (upliftLocationAndResult.UpliftFactor < upliftCriterion);
                double xCoordinateLastUpliftPoint = isUplift
                    ? upliftLocationAndResult.X
                    : currentSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X;
                var nonWaterRetaining2 = currentSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint2);
                if (nonWaterRetaining2 !=
                    null)
                {
                    xCoordinateLastUpliftPoint =
                        nonWaterRetaining2.X;
                }
                failureMechanismeParamatersMStab.MStabParameters.SlipCircleDefinition.XCoordinateLastUpliftPoint =
                    xCoordinateLastUpliftPoint;
            }

            failureMechanismeParamatersMStab.MStabParameters.ProjectFileName = mstabProjectFilename;
            failureMechanismeParamatersMStab.MStabParameters.SoilDatabaseName = scenario.Location.SoildatabaseName;

            if (!failureMechanismeParamatersMStab.IsComplete)
            {
                throw new Exception("Not all required data is available");
            }
            DamMStabAssembler assembler = new DamMStabAssembler();
            XDocument mstabXML = assembler.CreateDataTransferObject(failureMechanismeParamatersMStab);
            return mstabXML;
        }

        /// <summary>
        /// Determines which MStabResults of 2, contains the smallest safety factor
        /// </summary>
        /// <param name="mStabResults1"></param>
        /// <param name="mStabResults2"></param>
        /// <returns>MStab results with the minimum safety factor</returns>
        public static MStabResults MinMStabResults(MStabResults mStabResults1, MStabResults mStabResults2)
        {
            if (mStabResults2.zone1.safetyFactor < mStabResults1.zone1.safetyFactor)
            {
                return mStabResults2;
            }
            else
            {
                return mStabResults1;
            }
        }

        /// <summary>
        /// This is the place where we determine where the DGeoStability executable is located
        /// </summary>
        public static string MStabExePath
        {
            get
            {
                return Path.Combine(Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location),
                    @"DGeoStability.exe");
            }
        }

        public string StabilityProjectFilename
        {
            get { return stabilityProjectFilename; }
        }

        public List<LogMessage> ErrorMessages
        {
            get { return errorMessages; }
            set { errorMessages = value; }
        }

        // option of various computation kernel types
        public StabilityKernelType SelectedStabilityKernelType { get; set; }
        
        /// <summary>
        /// save stability project to dsx file
        /// </summary>
        /// <param name="fileName"></param>
        /// <param name="project"></param>
        public void SaveStabilityProjectToDsxFile(string fileName, MStabProject project)
        {
            DataEventPublisher.InvokeWithoutPublishingEvents(() =>
            {
                var xmlSerializer = new XmlSerializer();
                xmlSerializer.Serialize(project, fileName, AppDomain.CurrentDomain.BaseDirectory + "Mstab.xsl");
            });
        }

        private string SaveStabilityProjectFileToString(object stabilityProject)
        {
            var xmlSerializer = new XmlSerializer();
            return xmlSerializer.SerializeToString(stabilityProject);
        }

        /// <summary>
        /// Gets the location with lowest uplift factor.
        /// </summary>
        /// <param name="surfaceLine">The surface line.</param>
        /// <param name="soilProfile">The soil profile.</param>
        /// <param name="soilGeometry2DName">Name of the soil geometry2 D.</param>
        /// <param name="plLines">The pl lines.</param>
        /// <param name="location">The location.</param>
        /// <returns></returns>
        public UpliftLocationAndResult GetLocationWithLowestUpliftFactor(SurfaceLine2 surfaceLine,
            SoilProfile soilProfile, string soilGeometry2DName, PLLines plLines, Location location)
        {
            var profile1D = soilProfile as SoilProfile1D;
            var profile2D = soilProfile as SoilProfile2D;
            UpliftLocationDeterminator upliftLocationDeterminator = new UpliftLocationDeterminator()
            {
                SurfaceLine = surfaceLine,
                SoilProfile = profile1D,
                SoilProfile2D = profile2D,
                SoilGeometry2DName = soilGeometry2DName,
                SoilBaseDB = location.SoilbaseDB,
                SoilList = location.SoilList,
                DikeEmbankmentMaterial = location.GetDikeEmbankmentSoil(),
                PLLines = plLines,
                XSoilGeometry2DOrigin = location.XSoilGeometry2DOrigin
            };
            return upliftLocationDeterminator.GetLocationAtWithLowestUpliftFactor();
        }

        public void Dispose()
        {
            if (dotNetMstabProjectResults != null)
            {
                dotNetMstabProjectResults.Dispose();
            }
        }
    }
}