//-----------------------------------------------------------------------
// <copyright file="DamFailureMechanismeCalculator.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>16-12-2009</date>
// <summary>Perform single calculation which is specified in DamFailureMechanismeCalculationSpecification</summary>
//-----------------------------------------------------------------------

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Linq;
using System.Reflection;
using System.Text.RegularExpressions;
using System.Xml.Linq;
using Deltares.FlowSlide.Data;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Mathematics;
using Deltares.Piping.Data;
using Deltares.Probabilistic;
using Deltares.Standard;
using Deltares.Standard.IO.Xml;
using Deltares.Standard.Language;
using Deltares.Standard.Logging;
using Deltares.Standard.Validation;

namespace Deltares.Dam.Data
{
    public class MaximumRedesignIterationsReachedException : Exception
    {
        public MaximumRedesignIterationsReachedException() : base("Maximum number of iterations in redesign reached.")
        {
        }
    }

    /// <summary>
    /// 
    /// </summary>
    public class DamFailureMechanismeCalculator
    {
        private const int PipingRedesingMaxIterations = 400;
        private const double defaultMaxFractionOfDikeHeightForShoulderHeight = 0.67;
        private const double CMinimumShoulderElevation = 0.5;
        private const double CMinimumShoulderExtraElevation = 0.05;
        private const double CMinimumShoulderWidth = 2;
        private const double CMinimumShoulderExtraWidth = 1;
        private const double CToleranceShoulderChanges = 0.001;

        private readonly DamFailureMechanismeCalculationSpecification damFailureMechanismeCalculationSpecification;
        private readonly string mstabProgramPath;
        private readonly ProgramType programType;
        private readonly string slopeWProgramPath;
        private AnalysisType analysisType;
        private List<LogMessage> errorMessages = new List<LogMessage>();
        private int progressSofar;

        public DamFailureMechanismeCalculator(ProgramType programType, DamFailureMechanismeCalculationSpecification damFailureMechanismeCalculationSpecification,
                                              string mstabProgramPath, string slopeWProgramPath, string mapForSoilGeometries2D, ProbabilisticType probabilisticType)
        {
            this.programType = programType;
            this.damFailureMechanismeCalculationSpecification = damFailureMechanismeCalculationSpecification;
            this.mstabProgramPath = mstabProgramPath;
            this.slopeWProgramPath = slopeWProgramPath;
            MapForSoilGeometries2D = mapForSoilGeometries2D;
            ProbabilisticType = probabilisticType;
        }

        public string CalculationBaseDirectory { get; set; }
        public string MapForSoilGeometries2D { get; set; }
        public ProbabilisticType ProbabilisticType { get; set; }
        public NonWaterRetainingObject NonWaterRetainingObject { get; set; }

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

        /// <summary>
        /// Main calculation call
        /// </summary>
        /// <param name="analysisTypeGiven"></param>
        /// <param name="scenario"></param>
        public void Calculate(AnalysisType analysisTypeGiven, Scenario scenario)
        {
            Debug.WriteLine("Location '{0}', scenario '{1}' DamFailureMechanismeCalculator.Calculation", scenario.Location.Name, scenario.LocationScenarioID);
            analysisType = analysisTypeGiven;

            switch (analysisType)
            {
                case AnalysisType.NoAdaption:
                    CalculateForScenario(scenario);
                    break;

                case AnalysisType.AdaptGeometry:
                    RedesignSurfaceLinesForAllScenarios(scenario);
                    break;
                case AnalysisType.AdaptNWO:
                    ImplementNWOInSurfaceLinesForAllScenarios(scenario);
                    break;
            }
        }

        /// <summary>
        /// Selects the probabilities for specific failure mechanism.
        /// </summary>
        /// <param name="failureMechanismSystemType">Type of the failure mechanism system.</param>
        /// <param name="soilProfileProbabilities">The soil profile probabilities.</param>
        /// <returns></returns>
        public static List<SoilGeometryProbability> SelectProbabilitiesForFailureMechanism(FailureMechanismSystemType failureMechanismSystemType, IList<SoilGeometryProbability> soilProfileProbabilities)
        {
            switch (failureMechanismSystemType)
            {
                case FailureMechanismSystemType.StabilityInside:
                case FailureMechanismSystemType.StabilityOutside:
                case FailureMechanismSystemType.HorizontalBalance:
                    return SelectStabilityProbabilities(soilProfileProbabilities);
                case FailureMechanismSystemType.Piping:
                    return SelectPipingProbabilities(soilProfileProbabilities);
                case FailureMechanismSystemType.FlowSlide:
                    return SelectFlowSlideProbabilities(soilProfileProbabilities);
            }
            return null;
        }

        /// <summary>
        /// Redesigns the height of the surface line.
        /// </summary>
        /// <param name="failureMechanismSystemType">Type of the failure mechanism system.</param>
        /// <param name="scenario">The scenario.</param>
        /// <param name="surfaceLine">The surface line.</param>
        /// <returns></returns>
        public static SurfaceLine2 RedesignSurfaceLineHeight(FailureMechanismSystemType failureMechanismSystemType, Scenario scenario, SurfaceLine2 surfaceLine)
        {
            double? dikeHeight = surfaceLine.GetDikeHeight();
            if (dikeHeight.HasValue)
            {
                if (scenario.DikeTableHeight > dikeHeight.Value)
                {
                    var surfaceLineHeightAdapter = new SurfaceLineHeightAdapter(surfaceLine, scenario.Location);
                    SurfaceLine2 adaptedSurfaceLine = surfaceLineHeightAdapter.ConstructNewSurfaceLine(scenario.DikeTableHeight ?? surfaceLine.GetDefaultDikeTableHeight() ?? 0);

                    if (adaptedSurfaceLine != null)
                    {
                        var validationError = adaptedSurfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                        if (validationError != null)
                        {
                            adaptedSurfaceLine.Dispose();
                            throw new SurfaceLineException(validationError.Text);
                        }

                        surfaceLine = adaptedSurfaceLine;
                        foreach (var soilProfileProbability in scenario.Location.Segment.SoilProfileProbabilities)
                        {
                            if (soilProfileProbability.SegmentFailureMechanismType == failureMechanismSystemType || soilProfileProbability.SegmentFailureMechanismType == null)
                            {
                                scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, surfaceLine);
                            }
                        }
                    }
                }
            }
            return surfaceLine;
        }

        /// <summary>
        /// Cleanup directory where the piping project files are created
        /// </summary>
        public static void ClearPipingCalculationBaseDirectory()
        {
            string pipingBaseDirectory = GetPipingCalculationBaseDirectory();
            if (Directory.Exists(pipingBaseDirectory))
            {
                Directory.Delete(pipingBaseDirectory, true);
            }
        }

        private void CalculateForScenario(Scenario scenario)
        {
            switch (damFailureMechanismeCalculationSpecification.FailureMechanismSystemType)
            {
                case FailureMechanismSystemType.StabilityInside:
                case FailureMechanismSystemType.StabilityOutside:
                    CalculateStabilityForScenario(scenario, 0);
                    break;
                case FailureMechanismSystemType.Piping:
                    CalculatePipingForScenario(scenario);
                    break;
                case FailureMechanismSystemType.FlowSlide:
                    CalculateFlowSlideForScenario(scenario);
                    break;
            }
        }

        /// <summary>
        /// Checks the and adapt betas water levels when needed.
        /// </summary>
        /// <param name="betas">The betas.</param>
        /// <param name="waterLevels">The water levels.</param>
        /// <param name="scenario">The scenario.</param>
        /// <param name="newSurfaceLine">The new surface line.</param>
        /// <param name="pipingCalculator">The piping calculator.</param>
        /// <param name="soilProfileProbability">The soil profile probability.</param>
        /// <returns></returns>
        private bool CheckAndAdaptBetasWaterLevelsWhenNeeded(double[] betas, double[] waterLevels, Scenario scenario, SurfaceLine2 newSurfaceLine, PipingCalculator pipingCalculator, SoilGeometryProbability soilProfileProbability)
        {
            int nrValidItems = 0;
            for (int i = 0; i < 3; i++)
            {
                // If beta >= 10, then in fact no uplift occured and probability of failure is in fact 0. Such a value would be invalid for the design point calculation.
                if (betas[i] < 10)
                {
                    nrValidItems++;
                }
            }
            if ((nrValidItems == 0) || (nrValidItems == 3))
            {
                return (nrValidItems == 3);
            }
            else
            {
                if (nrValidItems == 1)
                {
                    // extra calculation(s) needed to determine a second valid beta (<10) to be able to perform a design point calculation. Search between toplevel 
                    // (beta[0]) and middle level a waterlevel at which a beta is valid and use that one (adapt the levels too of course). 
                    double? betaloc = null;
                    bool lReady = false;
                    double locLevel = (waterLevels[0] + waterLevels[1])*0.5;
                    var waterLevelProbabilistic = new ProbabilisticStruct(0.0, 0.01, (int) DistributionType.LogNormal, false);
                    waterLevelProbabilistic.Distribution = (int) DistributionType.LogNormal;
                    waterLevelProbabilistic.StandardDeviation = 0.01;
                    waterLevelProbabilistic.IsUsed = false;
                    while (!lReady)
                    {
                        waterLevelProbabilistic.Mean = locLevel;
                        betaloc = CalculatePipingReliabilityIndexForSurface(scenario, newSurfaceLine, pipingCalculator, soilProfileProbability,
                                                                            waterLevelProbabilistic);
                        lReady = betaloc.Value < 10;
                        if (lReady)
                        {
                            betas[1] = betaloc.Value;
                            waterLevels[1] = locLevel;
                        }
                        else
                        {
                            locLevel = locLevel + (waterLevels[0] + locLevel)*0.5;
                            lReady = Math.Abs(locLevel - waterLevels[0]) < 0.0001;
                        }
                    }

                    // now set nrValidItems to 2 to get the required third valid beta by axtrapolation
                    nrValidItems = 2;
                }
                if (nrValidItems == 2)
                {
                    // use the 2 valid betas to extrapolate a valid third beta value. 
                    betas[2] = Routines2D.LinInpolY(waterLevels[0], betas[0], waterLevels[1], betas[1], waterLevels[2]);
                }
                return true;
            }
        }

        private double? CalculatePipingFailureProbabilityAdvanced(PipingCalculator pipingCalculator, Scenario scenario,
                                                                  SoilGeometryProbability soilProfileProbability, SurfaceLine2 newSurfaceLine)
        {
            double[] waterLevels = scenario.DetermineProperWaterlevelsForProbabilisticAdvanced();
            var betas = new double[3];
            var waterLevelProbabilistic = new ProbabilisticStruct(0.0, 0.01, (int) DistributionType.LogNormal, false);
            for (int i = 0; i < 3; i++)
            {
                betas[i] = 1;
                waterLevelProbabilistic.Mean = waterLevels[i];
                string calculationName = String.Format("Prob{0}_Loc({1})_Sce({2})_Pro({3})_wl({4})",
                                                       PipingModelType.Sellmeijer.ToString(), scenario.Location.Name, scenario.LocationScenarioID, soilProfileProbability.SoilGeometryName, i);
                calculationName = Regex.Replace(calculationName, @"[\\\/:\*\?""'<>|.]", "_");
                pipingCalculator.FilenameCalculation = pipingCalculator.PipingCalculationDirectory + calculationName;

                double? betaloc = CalculatePipingReliabilityIndexForSurface(scenario, newSurfaceLine, pipingCalculator, soilProfileProbability, waterLevelProbabilistic);
                if (betaloc != null)
                {
                    betas[i] = betaloc.Value;
                }
            }

            // 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! 
            // After consulting Best, the next decision is made. MHW should be the riverlevel. And if the maxwaterlevel is smaller than MHW,then a designpoint
            // calculation is not longer required as the safety is determined by the maxwaterlevel. So just return the Beta[0] (= beta at maxwaterlevel) in 
            // that case.
            if (scenario.RiverLevel > scenario.MaxWaterLevel)
            {
                return Probabilistic.Probabilistic.NormalDistribution(-betas[0]);
            }
            else
            {
                if (CheckAndAdaptBetasWaterLevelsWhenNeeded(betas, waterLevels, scenario, newSurfaceLine, pipingCalculator, soilProfileProbability))
                {
                    var designPointWater = new DesignPointCalculation();
                    designPointWater.Betas = betas;
                    designPointWater.Waterlevels = waterLevels;
                    designPointWater.MHW = (double) (scenario.RiverLevel);
                    designPointWater.Decimate = (double) scenario.WaterHeightDecimeringsHoogte;
                    designPointWater.Exceed = DesignPointCalculation.ExceedingSet.twoThousend;
                    designPointWater.IsMaxLevelUsed = false;
                    designPointWater.MaxLevel = 0;
                    if (designPointWater.CalculateTheWaterDesignpoint())
                    {
                        return Probabilistic.Probabilistic.NormalDistribution(-designPointWater.Beta);
                    }
                    else
                    {
                        // result could not be determined.
                        throw new DamFailureMechanismeCalculatorException("Water design point calculation failed");
                    }
                }
                else
                {
                    // total situation is safe, so return 0.
                    return 0;
                }
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="newSurfaceLine"></param>
        /// <param name="pipingCalculator"></param>
        /// <param name="soilProfileProbability"></param>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        private double? CalculatePipingReliabilityIndexForSurface(Scenario scenario, SurfaceLine2 newSurfaceLine, PipingCalculator pipingCalculator, SoilGeometryProbability soilProfileProbability,
                                                                  ProbabilisticStruct waterLevel)
        {
            double? betaloc;
            if (newSurfaceLine != null)
            {
                betaloc = pipingCalculator.CalculateReliabilityIndex(scenario.Location, newSurfaceLine, soilProfileProbability.SoilProfile, waterLevel);
            }
            else
            {
                betaloc = pipingCalculator.CalculateReliabilityIndex(scenario.Location, scenario.GetMostRecentSurfaceLine(soilProfileProbability.SoilProfile,
                                                                                                                          soilProfileProbability.SoilGeometry2DName), soilProfileProbability.SoilProfile, waterLevel);
            }
            return betaloc;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="pipingCalculationDirectory"></param>
        /// <param name="scenario"></param>
        /// <param name="waterLevelProbabilistic"></param>
        private void CalculatePipingForScenario(Scenario scenario)
        {
            var waterLevelProbabilistic = new ProbabilisticStruct(0.0, 0.01, (int) DistributionType.LogNormal, false);
            string pipingCalculationDirectory = GetPipingCalculationBaseDirectory();
            try
            {
                var validSoilProfileProbabilities = SelectPipingProbabilities(scenario.Location.Segment.SoilProfileProbabilities);
                if (validSoilProfileProbabilities.Count == 0)
                {
                    throw new DamFailureMechanismeCalculatorException(String.Format("No piping profiles to calcutate for location {0}", scenario.Location.Name));
                }
                foreach (var soilProfileProbability in validSoilProfileProbabilities)
                {
                    try
                    {
                        var pipingResults = new PipingResults();
                        foreach (PipingModelType pipingCalculationType in Enum.GetValues(typeof(PipingModelType)))
                        {
                            var pipingProbabilisticParameters = new PipingProbabilisticParameters(scenario.Location.LayerHeightDistribution, scenario.Location.LayerHeightDeviation);
                            PipingCalculator pipingCalculator = pipingCalculatorFactory(scenario.Location, scenario.Location.CreateModelParametersForPLLines(),
                                                                                        pipingCalculationType,
                                                                                        scenario.GetUpliftCriterionPiping(scenario.Location.ModelFactors.UpliftCriterionPiping),
                                                                                        scenario.GetRequiredSafetyFactorPiping(scenario.Location.ModelFactors.RequiredSafetyFactorPiping),
                                                                                        pipingProbabilisticParameters);
                            pipingCalculator.PipingCalculationDirectory = pipingCalculationDirectory;
                            double? safetyFactor = null;
                            double? failureProbability = null;
                            double? failureProbabilityAdvanced = null;
                            double? pipingExitPointX = null;
                            double? upliftFactor = null;
                            double? heaveFactor = null;
                            try
                            {
                                switch (ProbabilisticType)
                                {
                                    case ProbabilisticType.Deterministic:
                                        safetyFactor = pipingCalculator.CalculatePipingFactor(
                                            scenario.Location,
                                            scenario.GetMostRecentSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName),
                                            soilProfileProbability.SoilProfile,
                                            scenario.RiverLevel);
                                        if (pipingCalculator.UpliftLocationAndResult != null)
                                        {
                                            pipingExitPointX = pipingCalculator.UpliftLocationAndResult.X;
                                            upliftFactor = pipingCalculator.UpliftLocationAndResult.UpliftFactor;
                                            heaveFactor = pipingCalculator.HeaveFactor;
                                        }
                                        break;
                                    case ProbabilisticType.Probabilistic:
                                        waterLevelProbabilistic.Mean = scenario.RiverLevel;

                                        string calculationName = String.Format("Prob{0}_Loc({1})_Sce({2})_Pro({3})",
                                                                               pipingCalculationType.ToString(), scenario.Location.Name, scenario.LocationScenarioID, soilProfileProbability.SoilGeometryName);
                                        calculationName = Regex.Replace(calculationName, @"[\\\/:\*\?""'<>|.]", "_");
                                        pipingCalculator.FilenameCalculation = pipingCalculator.PipingCalculationDirectory + calculationName;

                                        failureProbability =
                                            pipingCalculator.CalculatePipingFailureProbability(
                                                scenario.Location,
                                                scenario.GetMostRecentSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName),
                                                soilProfileProbability.SoilProfile,
                                                waterLevelProbabilistic);
                                        break;
                                    case ProbabilisticType.ProbabilisticFragility:
                                        failureProbabilityAdvanced = CalculatePipingFailureProbabilityAdvanced(pipingCalculator, scenario, soilProfileProbability, null);
                                        break;
                                }
                            }
                            catch (Exception exception)
                            {
                                throw new DamFailureMechanismeCalculatorException(exception.Message);
                            }

                            if (pipingCalculationType == PipingModelType.Sellmeijer)
                            {
                                scenario.SetSafetyFactorPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, safetyFactor);
                                if (ProbabilisticType == ProbabilisticType.ProbabilisticFragility)
                                {
                                    scenario.SetFailureProbabilityPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName,
                                                                         failureProbabilityAdvanced);
                                }
                                else
                                {
                                    scenario.SetFailureProbabilityPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName,
                                                                         failureProbability);
                                }
                            }
                            pipingResults.CalculationName = pipingCalculator.FilenameCalculation;
                            pipingResults.CalculationSubDir = @"PipingCalculation\\";
                            pipingResults.PipingExitPointX = pipingExitPointX;
                            pipingResults.UpliftFactor = upliftFactor;
                            pipingResults.HeaveFactor = heaveFactor;
                            switch (pipingCalculationType)
                            {
                                case PipingModelType.Bligh:
                                    pipingResults.BlighHCritical = pipingCalculator.HCritical;
                                    pipingResults.BlighPipingFactor = safetyFactor;
                                    break;
                                case PipingModelType.Sellmeijer:
                                    pipingResults.SellmeijerHCritical = pipingCalculator.HCritical;
                                    pipingResults.SellmeijerPipingFactor = safetyFactor;
                                    break;
                                case PipingModelType.Sellmeijer2Forces:
                                    pipingResults.Sellmeijer2ForcesHCritical = pipingCalculator.HCritical;
                                    pipingResults.Sellmeijer2ForcesPipingFactor = safetyFactor;
                                    break;
                                case PipingModelType.Sellmeijer4Forces:
                                    pipingResults.Sellmeijer4ForcesHCritical = pipingCalculator.HCritical;
                                    pipingResults.Sellmeijer4ForcesPipingFactor = safetyFactor;
                                    break;
                                case PipingModelType.Wti2017:
                                    pipingResults.Wti2017HCritical = pipingCalculator.HCritical;
                                    pipingResults.Wti2017PipingFactor = safetyFactor;
                                    break;
                            }
                        }
                        scenario.SetPipingResults(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, pipingResults);
                    }
                    catch (DamFailureMechanismeCalculatorException calculatorException)
                    {
                        string errorMessage = "FAIL: " + calculatorException.Message;
                        Exception innerException = calculatorException.InnerException;
                        while (innerException != null)
                        {
                            errorMessage = errorMessage + ";" + innerException.Message;
                            innerException = innerException.InnerException;
                        }
                        scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, errorMessage);

                        calculatorException.Scenario = scenario;
                        throw calculatorException;
                    }
                }
            }
            catch (Exception exception)
            {
                scenario.Errors.Add(exception.Message);
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    scenario.Errors.Add(innerException.Message);
                    innerException = innerException.InnerException;
                }
            }
        }

        /// <summary>
        /// Calculates failuremechanism flowslide for scenario.
        /// </summary>
        /// <param name="scenario">The scenario.</param>
        private void CalculateFlowSlideForScenario(Scenario scenario)
        {
            try
            {
                var validSoilProfileProbabilities = SelectFlowSlideProbabilities(scenario.Location.Segment.SoilProfileProbabilities);
                if (validSoilProfileProbabilities.Count == 0)
                {
                    throw new DamFailureMechanismeCalculatorException(String.Format("No profiles to calcutate for location {0}", scenario.Location.Name));
                }
                foreach (var soilProfileProbability in validSoilProfileProbabilities)
                {
                    try
                    {
                        if (soilProfileProbability.SoilProfile == null)
                        {
                            throw new DamFailureMechanismeCalculatorException(String.Format("Soilprofile is not 1D for location {0}",
                                                                                            scenario.Location.Name));
                        }
                        var flowSlideCalculations = new FlowSlideCalculations();
                        using (var flowSlideProject = CreateFlowSlideProject(scenario, soilProfileProbability))
                        {
                            string calculationName = StabilityCalculator.DetermineCalculationFilename(scenario.Location.Name, scenario.LocationScenarioID,
                                                                                                      soilProfileProbability.SoilGeometryName, 0);
                            WriteFlowSlideProjectToFile(flowSlideProject, calculationName);

                            flowSlideCalculations.RunCalculation(flowSlideProject);
                            double safetyFactor = flowSlideProject.ResultOverallFactor;
                            scenario.SetSafetyFactorFlowSlide(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, safetyFactor);
                        }
                    }
                    catch (DamFailureMechanismeCalculatorException calculatorException)
                    {
                        string errorMessage = "FAIL: " + calculatorException.Message;
                        Exception innerException = calculatorException.InnerException;
                        while (innerException != null)
                        {
                            errorMessage = errorMessage + ";" + innerException.Message;
                            innerException = innerException.InnerException;
                        }
                        scenario.SetResultMessage(soilProfileProbability.SoilProfile,
                                                  soilProfileProbability.SoilGeometry2DName, errorMessage);

                        calculatorException.Scenario = scenario;
                        throw calculatorException;
                    }
                }
            }
            catch (Exception exception)
            {
                scenario.Errors.Add(exception.Message);
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    scenario.Errors.Add(innerException.Message);
                    innerException = innerException.InnerException;
                }
            }
        }

        /// <summary>
        /// Creates the flowslide project.
        /// </summary>
        /// <param name="scenario">The scenario.</param>
        /// <param name="soilProfileProbability">The soil profile probability.</param>
        /// <returns></returns>
        private static FlowSlideProject CreateFlowSlideProject(Scenario scenario, SoilGeometryProbability soilProfileProbability)
        {
            var flowSlideProject = new FlowSlideProject();
            flowSlideProject.FlowSlideCalculationParameters.UseFlowSlideCalculation = UseFlowSlideCalculation.Global;
            flowSlideProject.SurfaceLine2 = scenario.Location.LocalXZSurfaceLine2;
            flowSlideProject.SoilProfile = soilProfileProbability.SoilProfile;
            flowSlideProject.Soils.Clear();
            var projectSoils = new List<Soil>();
            foreach (var layer in soilProfileProbability.SoilProfile.Layers)
            {
                if (!projectSoils.Exists(T => layer.Soil.Name.Equals(T.Name)))
                {
                    projectSoils.Add(layer.Soil);
                    flowSlideProject.Soils.Add(layer.Soil);
                }
            }
            flowSlideProject.FlowSlideCalculationParameters.CalculationParametersGeneral.Waterlevel = scenario.RiverLevel;
            return flowSlideProject;
        }

        /// <summary>
        /// Writes the flowslide project to file.
        /// </summary>
        /// <param name="flowSlideProject">The flow slide project.</param>
        /// <param name="calculationName">Name of the calculation.</param>
        private void WriteFlowSlideProjectToFile(FlowSlideProject flowSlideProject, string calculationName)
        {
            string flowslideCalculationDirectory = Path.Combine(CalculationBaseDirectory, @"Flowslide\");
            if (!Directory.Exists(flowslideCalculationDirectory))
            {
                Directory.CreateDirectory(flowslideCalculationDirectory);
            }
            string filenameExtension = ".fsx";
            string fileName = calculationName + filenameExtension;
            string flowslideProjectFilename = Path.Combine(flowslideCalculationDirectory, fileName);
            var xmlSerializer = new XmlSerializer();
            xmlSerializer.Serialize(flowSlideProject, flowslideProjectFilename);
        }

        /// <summary>
        /// Selects the soilProfileProbabilities for piping.
        /// </summary>
        /// <param name="soilProfileProbabilities">The soil profile probabilities.</param>
        /// <returns></returns>
        private static List<SoilGeometryProbability> SelectPipingProbabilities(IList<SoilGeometryProbability> soilProfileProbabilities)
        {
            var validSoilProfileProbabilities =
                new List<SoilGeometryProbability>(soilProfileProbabilities.Where(
                    p => ((p.SegmentFailureMechanismType == FailureMechanismSystemType.Piping) ||
                          (p.SegmentFailureMechanismType == null))));
            return validSoilProfileProbabilities;
        }

        /// <summary>
        /// Selects the soilProfileProbabilities for stability.
        /// </summary>
        /// <param name="soilProfileProbabilities">The soil profile probabilities.</param>
        /// <returns></returns>
        private static List<SoilGeometryProbability> SelectStabilityProbabilities(IList<SoilGeometryProbability> soilProfileProbabilities)
        {
            var validSoilProfileProbabilities = new List<SoilGeometryProbability>(soilProfileProbabilities.Where(
                p => ((p.SegmentFailureMechanismType == FailureMechanismSystemType.StabilityInside) ||
                      (p.SegmentFailureMechanismType == FailureMechanismSystemType.StabilityOutside) ||
                      (p.SegmentFailureMechanismType == null))));
            return validSoilProfileProbabilities;
        }

        /// <summary>
        /// Selects the flowslide probabilities.
        /// </summary>
        /// <param name="soilProfileProbabilities">The soil profile probabilities.</param>
        /// <returns></returns>
        private static List<SoilGeometryProbability> SelectFlowSlideProbabilities(IList<SoilGeometryProbability> soilProfileProbabilities)
        {
            // For now use the stability probabilities
            return SelectStabilityProbabilities(soilProfileProbabilities);
        }

        /// <summary>
        /// Consistency check to be performed before stability calculation
        /// </summary>
        /// <param name="scenarios"></param>
        private void ConsistencyCheckStability(Scenario scenario)
        {
            ConsistencyCheckStabilityPerScenario(scenario);
        }

        private void ConsistencyCheckStabilityPerScenario(Scenario scenario)
        {
            // If outward stability calculation then low waterlevel is required
            if (damFailureMechanismeCalculationSpecification.FailureMechanismeParamatersMStab.MStabParameters.GridPosition == MStabGridPosition.Left)
            {
                if (!scenario.RiverLevelLow.HasValue)
                {
                    throw new DamFailureMechanismeCalculatorException(String.Format(
                        "Location {0} scenario {1} has no low waterlevel defined, which is required for outward stability calculation (grid at left).",
                        scenario.Location.Name, scenario.LocationScenarioID));
                }
            }
        }

        /// <summary>
        /// Stability calculation
        /// </summary>
        /// <param name="scenario"></param>
        private void CalculateStabilityForScenario(Scenario scenario, int iter)
        {
            try
            {
                ConsistencyCheckStability(scenario);
                Location location = scenario.Location;
                var modelParametersForPLLines = new ModelParametersForPLLines
                {
                    PenetrationLength = location.PenetrationLength,
                    DampingFactorPL3 = location.DampingFactorPL3,
                    DampingFactorPL4 = location.DampingFactorPL4,
                    PLLineCreationMethod = location.PLLineCreationMethod
                };
                var requiredSafetyFactorStabilityInnerSlope = scenario.RequiredSafetyFactorStabilityInnerSlope ??
                                                              location.ModelFactors.RequiredSafetyFactorStabilityInnerSlope;
                using (var stabilityCalculator =
                    new StabilityCalculator(
                        damFailureMechanismeCalculationSpecification.FailureMechanismeParamatersMStab,
                        programType,
                        modelParametersForPLLines,
                        location.TrafficLoad,
                        location.MinimalCircleDepth,
                        requiredSafetyFactorStabilityInnerSlope.Value,
                        mstabProgramPath, slopeWProgramPath,
                        location.GaugePLLines,
                        location.Gauges,
                        location.SoilbaseDB,
                        location.SoilList,
                        ProbabilisticType)
                    {
                        SelectedStabilityKernelType = damFailureMechanismeCalculationSpecification.StabilityKernelType
                    })
                {
                    //stabilityCalculator.CalculationBaseDirectory = this.CalculationBaseDirectory;
                    var validSoilProfileProbabilities = SelectStabilityProbabilities(scenario.Location.Segment.SoilProfileProbabilities);
                    if (validSoilProfileProbabilities.Count == 0)
                    {
                        throw new DamFailureMechanismeCalculatorException(String.Format("No stability profiles to calcutate for location {0}", scenario.Location.Name));
                    }
                    int errorCount = 0;
                    var scenarioErrorMessages = new List<String>();
                    foreach (var soilProfileProbability in validSoilProfileProbabilities)
                    {
                        try
                        {
                            UpliftSituation? upliftSituation = scenario.GetStabilityUpliftSituation(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName);
                            // The calculation for UpliftVan will only will done if there is Uplift 
                            if (IsCalculationRequired(upliftSituation))
                            {
                                stabilityCalculator.Calculate(scenario, soilProfileProbability.SoilProfile, GetFullSoilGeometry2DName(soilProfileProbability.SoilGeometry2DName), iter);
                                MStabResults? mStabResults = scenario.GetMStabResults(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName);
                                if ((mStabResults != null) && (AnalysisType.AdaptNWO == analysisType))
                                {
                                    // Store the NWO results
                                    stabilityCalculator.NWOPhreaticAdaption = NonWaterRetainingObject.PhreaticAdaption;
                                    var results = new NonWaterRetainingObjectResults();
                                    results.NwoId = NonWaterRetainingObject.NwoId;
                                    var coordinateSystemConverter = new CoordinateSystemConverter();
                                    coordinateSystemConverter.DefineGlobalXYZBasedOnLine(scenario.Location.SurfaceLine2.Geometry);
                                    results.AdaptedSurfaceLine = scenario.Location.LocalXZSurfaceLine2.FullDeepClone();
                                    coordinateSystemConverter.ConvertLocalXZToGlobalXYZ(results.AdaptedSurfaceLine.Geometry);
                                    results.LocationXrdStart = results.AdaptedSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X;
                                    results.LocationYrdStart = results.AdaptedSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).Y;
                                    results.LocationZrdStart = results.AdaptedSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).Z;
                                    results.LocationXrdEnd = results.AdaptedSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).X;
                                    results.LocationYrdEnd = results.AdaptedSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).Y;
                                    results.LocationZrdEnd = results.AdaptedSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).Z;

                                    results.MStabResults = mStabResults.Value;
                                    results.SoilProfileProbability = soilProfileProbability;
                                    scenario.NwoResults.Add(results);
                                }
                            }
                            // assign original surfaceline to redesigned surfaceline (not for NWODesign as that would provide wrong surface line) 
                            if (analysisType != AnalysisType.AdaptNWO)
                            {
                                scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile,
                                                                  soilProfileProbability.SoilGeometry2DName,
                                                                  scenario.Location.SurfaceLine2);
                            }
                        }
                        catch (DamFailureMechanismeCalculatorException calculatorException)
                        {
                            string errorMessage = "FAIL: " + calculatorException.Message;
                            Exception innerException = calculatorException.InnerException;
                            while (innerException != null)
                            {
                                errorMessage = errorMessage + ";" + innerException.Message;
                                innerException = innerException.InnerException;
                            }
                            scenario.SetResultMessage(soilProfileProbability.SoilProfile,
                                                      soilProfileProbability.SoilGeometry2DName, errorMessage);

                            calculatorException.Scenario = scenario;
                            errorCount++;
                            scenarioErrorMessages.Add(errorMessage);
                        }
                    }
                    if (errorCount > 0)
                    {
                        string errorMessage = String.Format("{0} calculation error(s) in location {1} scenario {2}",
                                                            errorCount, location.Name, scenario.LocationScenarioID);
                        foreach (var scenarioErrorMessage in scenarioErrorMessages)
                        {
                            errorMessage = String.Format("{0}; {1}", errorMessage, scenarioErrorMessage);
                        }
                        var calculatorException = new DamFailureMechanismeCalculatorException(errorMessage)
                        {
                            Scenario = scenario
                        };
                        throw calculatorException;
                    }
                }
            }
            catch (Exception exception)
            {
                scenario.Errors.Add(exception.Message);
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    scenario.Errors.Add(innerException.Message);
                    innerException = innerException.InnerException;
                }
                throw;
            }
        }

        /// <summary>
        /// Check if calculation should be done
        /// (In case of UpliftVan and no uplift this is not the case)
        /// </summary>
        /// <param name="upliftSituation">The uplift situation.</param>
        /// <returns>
        ///   <c>true</c> if [is calculation required] ; otherwise, <c>false</c>.
        /// </returns>
        private bool IsCalculationRequired(UpliftSituation? upliftSituation)
        {
            bool isUplift = upliftSituation.HasValue && upliftSituation.Value.IsUplift;
            bool isSkipCalculation = (!isUplift &&
                                      damFailureMechanismeCalculationSpecification.FailureMechanismeParamatersMStab.IsStabilityCheckOnUplift &&
                                      (damFailureMechanismeCalculationSpecification.FailureMechanismSystemType == FailureMechanismSystemType.StabilityInside) &&
                                      (damFailureMechanismeCalculationSpecification.FailureMechanismeParamatersMStab.MStabParameters.Model == MStabModelType.UpliftVan));
            return !isSkipCalculation;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenarios"></param>
        private void ImplementNWOInSurfaceLinesForAllScenarios(Scenario scenario)
        {
            try
            {
                ImplementNWOInSurfaceLine(scenario);
            }
            catch (Exception exception)
            {
                scenario.Errors.Add(exception.Message);
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    scenario.Errors.Add(innerException.Message);
                    innerException = innerException.InnerException;
                }
            }
        }

        private void CalculateScenario(Scenario scenario, int iter)
        {
            switch (damFailureMechanismeCalculationSpecification.FailureMechanismSystemType)
            {
                case FailureMechanismSystemType.StabilityInside:
                case FailureMechanismSystemType.StabilityOutside:
                    ConsistencyCheckStabilityPerScenario(scenario);
                    CalculateStabilityForScenario(scenario, iter);
                    break;

                case FailureMechanismSystemType.Piping:
                    CalculatePipingForScenario(scenario);
                    break;
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenario"></param>
        /// <returns></returns>
        private double? ImplementNWOInSurfaceLine(Scenario scenario)
        {
            double? safetyFactor = null;

            if (scenario.Location != null)
            {
                SurfaceLine2 originalLocalXzSurfaceLine = scenario.Location.LocalXZSurfaceLine2.FullDeepClone();
                try
                {
                    // fit the non water retaining object in the current surfaceline. If it does (no longer) fit, surfaceline is returned as null and 
                    // the calculation is finished. The phreaticline may have to be adapted to the new surfaceline.

                    var nwoInSurfacleLine = new NonWaterRetainingObjectInSurfaceLine();
                    nwoInSurfacleLine.NonWaterRetainingObject = NonWaterRetainingObject;
                    nwoInSurfacleLine.SurfaceLine = scenario.Location.LocalXZSurfaceLine2;
                    nwoInSurfacleLine.GridPosition = damFailureMechanismeCalculationSpecification.FailureMechanismeParamatersMStab.MStabParameters.GridPosition;
                    nwoInSurfacleLine.StepSizeX = NonWaterRetainingObject.StepSizeX;
                    double fitPositionX = nwoInSurfacleLine.DetermineStartLocationForNonWaterRetainingObject(true);
                    if (nwoInSurfacleLine.DoesPositionXFitInSurfaceLine(fitPositionX))
                    {
                        var iter = 0;
                        bool doneFittingNWO = false;
                        while (!doneFittingNWO)
                        {
                            // try to fit NWO in surfaceline. If successfull, calculate it, else this is finished.
                            SurfaceLine2 localXzSurfaceLine = nwoInSurfacleLine.FitNonWaterRetainingObjectInSurfaceLine(fitPositionX);
                            if (localXzSurfaceLine != null)
                            {
                                var validationError = localXzSurfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                                if (validationError != null)
                                {
                                    localXzSurfaceLine.Dispose();
                                    throw new SurfaceLineException(validationError.Text);
                                }

                                scenario.Location.LocalXZSurfaceLine2.Dispose();
                                scenario.Location.LocalXZSurfaceLine2 = localXzSurfaceLine;
                                // calculate with this surfaceline and store the results.
                                CalculateScenario(scenario, iter);
                                // Get new fitPositionX (based on stepsize X and gridposition)
                                fitPositionX = nwoInSurfacleLine.DetermineNewFitPostionX(fitPositionX, true);
                                iter++;
                            }
                            else
                            {
                                doneFittingNWO = true;
                            }
                        }

                        if (iter > 0)
                        {
                            // restore the original surface line
                            // Note: 'originalLocalXzSurfaceLine' is only used if now exception have occurred and
                            scenario.Location.LocalXZSurfaceLine2 = originalLocalXzSurfaceLine;
                        }
                        else
                        {
                            throw new DamFailureMechanismeCalculatorException(
                                "NonWaterRetainingObject does not fit within Surfaceline.");
                        }
                    }
                    else
                    {
                        throw new DamFailureMechanismeCalculatorException(
                            "NonWaterRetainingObject does not fit within Surfaceline.");
                    }
                }
                catch (Exception e)
                {
                    originalLocalXzSurfaceLine.Dispose(); // Clone has become an orphan, so dispose and discard

                    // Add scenario if having a failure mechanism calculation exception:
                    var calculatorException = e as DamFailureMechanismeCalculatorException;
                    if (calculatorException != null)
                    {
                        calculatorException.Scenario = scenario;
                    }
                    throw;
                }
            }

            return safetyFactor;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenario"></param>
        private void RedesignSurfaceLinesForAllScenarios(Scenario scenario)
        {
            try
            {
                RedesignSurfaceLine(scenario);
            }
            catch (Exception exception)
            {
                scenario.Errors.Add(exception.Message);
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    scenario.Errors.Add(innerException.Message);
                    innerException = innerException.InnerException;
                }
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="scenario"></param>
        /// <returns></returns>
        private void RedesignSurfaceLine(Scenario scenario)
        {
            try
            {
                if (scenario.Location != null)
                {
                    SurfaceLine2 surfaceLine = scenario.Location.LocalXZSurfaceLine2.FullDeepClone(); //TODO: I really have no clue this object should be in a using clause or not... :(
                    if (scenario.Location.RedesignDikeHeight)
                    {
                        // Dike height adaptation

                        var redesignedSurfaceLine = RedesignSurfaceLineHeight(damFailureMechanismeCalculationSpecification.FailureMechanismSystemType,
                                                                              scenario, surfaceLine);
                        var validationError = redesignedSurfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                        if (validationError != null)
                        {
                            redesignedSurfaceLine.Dispose();
                            throw new SurfaceLineException(validationError.Text);
                        }

                        // No error, so replace clone with redesigned surfaceline
                        surfaceLine.Dispose();
                        surfaceLine = redesignedSurfaceLine;
                    }

                    if (scenario.Location.RedesignDikeShoulder)
                    {
                        // Shoulder and slope adaption
                        double? safetyFactor = null;
                        switch (damFailureMechanismeCalculationSpecification.FailureMechanismSystemType)
                        {
                            case FailureMechanismSystemType.StabilityInside:
                            case FailureMechanismSystemType.StabilityOutside:
                                safetyFactor = scenario.ModelFactors.RequiredSafetyFactorStabilityInnerSlope;
                                RedesignSurfaceLineStabilityInside(scenario, safetyFactor);
                                break;

                            case FailureMechanismSystemType.Piping:
                                safetyFactor = scenario.ModelFactors.RequiredSafetyFactorPiping;
                                RedesignSurfaceLinePipingShoulder(damFailureMechanismeCalculationSpecification.PipingModelType,
                                                                  scenario, ref safetyFactor, ref surfaceLine);
                                break;
                            case FailureMechanismSystemType.FlowSlide:
                                RedesignSurfaceLineFlowSlide(scenario, safetyFactor, ref surfaceLine);
                                break;
                        }
                    }
                }
            }
            catch (DamFailureMechanismeCalculatorException calculatorException)
            {
                calculatorException.Scenario = scenario;
                throw;
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="failureMechanismType"></param>
        /// <param name="scenario"></param>
        /// <param name="safetyFactor"></param>
        /// <param name="surfaceLine"></param>
        private void RedesignSurfaceLinePipingShoulder(PipingModelType pipingModelType, Scenario scenario, ref double? safetyFactor, ref SurfaceLine2 surfaceLine)
        {
            var modelParametersForPLLines = new ModelParametersForPLLines
            {
                PenetrationLength = scenario.Location.PenetrationLength,
                DampingFactorPL3 = scenario.Location.DampingFactorPL3,
                DampingFactorPL4 = scenario.Location.DampingFactorPL4,
                PLLineCreationMethod = scenario.Location.PLLineCreationMethod
            };

            var pipingProbabilisticParameters = new PipingProbabilisticParameters(scenario.Location.LayerHeightDistribution,
                                                                                  scenario.Location.LayerHeightDeviation);
            var pipingUpliftCriterion = scenario.GetUpliftCriterionPiping(scenario.Location.ModelFactors.UpliftCriterionPiping);
            if (pipingUpliftCriterion <= 0)
            {
                throw new DamFailureMechanismeCalculatorException(String.Format("Invalid piping uplift criterion {0} for location {1} scenario {2}",
                                                                                pipingUpliftCriterion, scenario.Location.Name, scenario.LocationScenarioID));
            }

            var requiredSafetyFactorPiping = 0.0;
            if (scenario.RequiredSafetyFactorPiping != null)
            {
                requiredSafetyFactorPiping = scenario.RequiredSafetyFactorPiping.Value;
            }
            if (safetyFactor != null)
            {
                requiredSafetyFactorPiping = safetyFactor.Value;
            }
            if (requiredSafetyFactorPiping <= 0)
            {
                throw new DamFailureMechanismeCalculatorException(String.Format("Invalid required safetyfactor piping {0} for location {1} scenario {2}",
                                                                                requiredSafetyFactorPiping, scenario.Location.Name, scenario.LocationScenarioID));
            }
            PipingCalculator pipingCalculator = pipingCalculatorFactory(scenario.Location, modelParametersForPLLines, pipingModelType,
                                                                        pipingUpliftCriterion,
                                                                        requiredSafetyFactorPiping,
                                                                        pipingProbabilisticParameters);

            var validSoilProfileProbabilities = SelectPipingProbabilities(scenario.Location.Segment.SoilProfileProbabilities);
            if (validSoilProfileProbabilities.Count == 0)
            {
                throw new DamFailureMechanismeCalculatorException(String.Format("No piping profiles to calcutate for location {0}", scenario.Location.Name));
            }
            foreach (var soilProfileProbability in validSoilProfileProbabilities)
            {
                try
                {
                    SurfaceLine2 newSurfaceLine = surfaceLine.FullDeepClone();
                    switch (ProbabilisticType)
                    {
                        case ProbabilisticType.Deterministic:
                        {
                            var newCandidate = DetermineNewSafeSurfaceLinePipingDeterministic(scenario, soilProfileProbability, pipingCalculator, newSurfaceLine, ref safetyFactor);
                            if (!ReferenceEquals(newSurfaceLine, newCandidate))
                            {
                                newSurfaceLine.Dispose();
                                newSurfaceLine = newCandidate;
                            }

                            scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, newSurfaceLine);
                            scenario.SetSafetyFactorPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, safetyFactor);
                            PipingResults pipingResults = GetPipingResults(safetyFactor, pipingModelType, pipingCalculator);
                            scenario.SetPipingResults(soilProfileProbability.SoilProfile,
                                                      soilProfileProbability.SoilGeometry2DName, pipingResults);
                            break;
                        }

                        case ProbabilisticType.Probabilistic:
                        case ProbabilisticType.ProbabilisticFragility:
                        {
                            //Todo Adapt the determination of the piping design to the way it should be (See Detrministic and functional design)
                            double[] waterLevels = scenario.DetermineProperWaterlevelsForProbabilisticAdvanced();
                            double waterLevel = scenario.RiverLevel;

                            // for test of deterministic as initialisation for probabilistic, use next line.
                            // newSurfaceLine = GetNewSafeSurfaceLine(scenario, soilProfileProbability, pipingCalculator, newSurfaceLine, ref requiredSafetyFactor);

                            bool upLiftOccured = false;
                            if (ProbabilisticType == ProbabilisticType.ProbabilisticFragility)
                            {
                                waterLevel = waterLevels[1];
                            }
                            // start by checking the uplift first (if so, a newShoulderHeight can be determined). If uplift is not an issue then 
                            // probability of failure = 0. If uplift is an issue, add a berm until uplift is no longer an issue.
                            double? newShoulderHeight = pipingCalculator.CalculateDesignShoulderHeight(scenario.Location, newSurfaceLine,
                                                                                                       soilProfileProbability.SoilProfile, waterLevel);

                            double currentShoulderHeight = 0;
                            double currentShoulderWidth = 0;
                            double newShoulderWidth = 0;
                            double slopeFactor = 0;
                            if (newShoulderHeight != null)
                            {
                                // Debug.Write(String.Format("New Shoulder Height {0:0.000}\n", newShoulderHeight));

                                // if a new shoulderheight could be found, uplift is an issue. The new found shoulderheight might however not be the correct one
                                // as with this height we insert or adjust a berm. This new/adjusted berm has to be checked until uplift is no longer an issue.
                                upLiftOccured = true;

                                // Start by adding the new shoulder. Note that there might be an existing shoulder.
                                var shoulderTopInside = newSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderTopInside);
                                var dikeToeAtPolder = newSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
                                if (newSurfaceLine.HasShoulderInside())
                                {
                                    // there was a shoulder which now has to be extended. 
                                    currentShoulderWidth = shoulderTopInside.X -
                                                           newSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside).X;
                                    currentShoulderHeight = shoulderTopInside.Z - dikeToeAtPolder.Z;
                                    // shoulder must grow in length, it should be passed the uplift location
                                    newShoulderWidth = currentShoulderWidth + pipingCalculator.UpliftLocationAndResult.X - shoulderTopInside.X;
                                    // in height it should grow with the calculated amount.
                                    newShoulderHeight = newShoulderHeight.Value + currentShoulderHeight;
                                }
                                else
                                {
                                    // there was no shoulder, so create one
                                    var dikeTopAtPolder = newSurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
                                    slopeFactor = (dikeToeAtPolder.X - dikeTopAtPolder.X)/(dikeTopAtPolder.Z - dikeToeAtPolder.Z);
                                    newShoulderWidth = 0.5 + newShoulderHeight.Value*slopeFactor;
                                }
                                // Make sure new berm is not too small
                                newShoulderHeight = Math.Max(newShoulderHeight.Value, CMinimumShoulderElevation);
                                newShoulderWidth = Math.Max(newShoulderWidth, CMinimumShoulderWidth);
                                currentShoulderWidth = newShoulderWidth;
                                currentShoulderHeight = newShoulderHeight.Value;

                                bool finished = false;
                                double maxShoulderLevel = CalculateMaximumShoulderLevel(newSurfaceLine, 1.0); // no limit to height of shoulder
                                while ((newShoulderHeight != null) && (!finished))
                                {
                                    var surfaceLineShoulderAdapter = new SurfaceLineShoulderAdapter(newSurfaceLine, scenario.Location);
                                    surfaceLineShoulderAdapter.MaxShoulderLevel = maxShoulderLevel;
                                    newSurfaceLine = surfaceLineShoulderAdapter.ConstructNewSurfaceLine(newShoulderWidth, newShoulderHeight.Value, true);
                                    newShoulderHeight = pipingCalculator.CalculateDesignShoulderHeight(scenario.Location, newSurfaceLine,
                                                                                                       soilProfileProbability.SoilProfile, waterLevel);
                                    finished = (pipingCalculator.UpliftLocationAndResult.X > shoulderTopInside.X);
                                    if (newShoulderHeight != null && (!finished))
                                    {
                                        // make sure new height is at least cMinimumShoulderExtraElevation higher
                                        if (newShoulderHeight.Value < CMinimumShoulderExtraElevation)
                                        {
                                            newShoulderHeight = currentShoulderHeight + CMinimumShoulderExtraElevation;
                                        }
                                        else
                                        {
                                            newShoulderHeight = currentShoulderHeight + newShoulderHeight.Value;
                                        }
                                        // shoulder must grow in length, it should be passed the uplift location and at least the minimum grow value
                                        newShoulderWidth = currentShoulderWidth + CMinimumShoulderExtraWidth;
                                        currentShoulderWidth = currentShoulderWidth + pipingCalculator.UpliftLocationAndResult.X -
                                                               shoulderTopInside.X;
                                        newShoulderWidth = Math.Max(currentShoulderWidth, newShoulderWidth);

                                        currentShoulderWidth = newShoulderWidth;
                                        currentShoulderHeight = newShoulderHeight.Value;
                                    }
                                }
                            }

                            // if no shoulder needed (i.e. when no uplift occurs) then failure probability is always cDefaultMinReturnValue = 0.0
                            double? failureProbability = PipingCalculator.cDefaultMinReturnValue;
                            if (upLiftOccured)
                            {
                                // Uplift occured, so calculate failureProbability
                                if (ProbabilisticType == ProbabilisticType.Probabilistic)
                                {
                                    var waterLevelProbabilistic = new ProbabilisticStruct(0.0, 0.01, (int) DistributionType.LogNormal, false);
                                    waterLevelProbabilistic.Mean = scenario.RiverLevel;
                                    failureProbability = pipingCalculator.CalculatePipingFailureProbability(scenario.Location, newSurfaceLine,
                                                                                                            soilProfileProbability.SoilProfile, waterLevelProbabilistic);
                                }
                                else
                                {
                                    failureProbability = CalculatePipingFailureProbabilityAdvanced(pipingCalculator, scenario, soilProfileProbability,
                                                                                                   newSurfaceLine);
                                }

                                if (failureProbability != null)
                                {
                                    int iterationIndex = 1;
                                    bool isRedesignRequired = (failureProbability.Value > scenario.ModelFactors.RequiredProbabilityOfFailurePiping);

                                    double maxShoulderLevel = CalculateMaximumShoulderLevel(newSurfaceLine, 1.0); // no limit to height of shoulder
                                    while (isRedesignRequired && newSurfaceLine != null)
                                    {
                                        // Debug.Write(String.Format("Iteration {0} failureProbability {1:0.0000000000}, currentShoulderWidth {2:0.000}, currentShoulderHeight {3:0.000}\n", iterationIndex, failureProbability, currentShoulderWidth, currentShoulderHeight));
                                        iterationIndex++;
                                        currentShoulderWidth = currentShoulderWidth + 1;
                                        // Due to different calculation types, raise the shoulder too next to making it wider.
                                        currentShoulderHeight = currentShoulderHeight + CMinimumShoulderExtraElevation;
                                        var surfaceLineShoulderAdapter = new SurfaceLineShoulderAdapter(newSurfaceLine, scenario.Location);
                                        surfaceLineShoulderAdapter.MaxShoulderLevel = maxShoulderLevel;
                                        newSurfaceLine = surfaceLineShoulderAdapter.ConstructNewSurfaceLine(currentShoulderWidth, currentShoulderHeight, true);
                                        if (ProbabilisticType == ProbabilisticType.Probabilistic)
                                        {
                                            var waterLevelProbabilistic = new ProbabilisticStruct(0.0, 0.01, (int) DistributionType.LogNormal, false);
                                            waterLevelProbabilistic.Mean = scenario.RiverLevel;
                                            failureProbability = pipingCalculator.CalculatePipingFailureProbability(
                                                scenario.Location, newSurfaceLine,
                                                soilProfileProbability.SoilProfile, waterLevelProbabilistic);
                                        }
                                        else
                                        {
                                            failureProbability = CalculatePipingFailureProbabilityAdvanced(pipingCalculator, scenario,
                                                                                                           soilProfileProbability, newSurfaceLine);
                                        }

                                        if (failureProbability != null)
                                        {
                                            isRedesignRequired = (failureProbability.Value > scenario.ModelFactors.RequiredProbabilityOfFailurePiping);
                                        }
                                        else
                                        {
                                            isRedesignRequired = true;
                                        }
                                        if (iterationIndex >= PipingRedesingMaxIterations)
                                        {
                                            // #BKA: vraag voor Ray set the results in case of failure too to provide feedback. (of niet Tom???) Voorlopig even niet. Ray was net naar huis dus nog navragen.
                                            //  scenario.SetFailureProbabilityPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName,
                                            //    failureProbability);
                                            // scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, 
                                            //   newSurfaceLine);
                                            throw new MaximumRedesignIterationsReachedException();
                                        }
                                    }
                                }
                            }
                            if (failureProbability != null)
                            {
                                scenario.SetFailureProbabilityPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName,
                                                                     failureProbability);
                            }
                            else
                            {
                                scenario.SetFailureProbabilityPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName,
                                                                     999);
                            }
                            scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, newSurfaceLine);
                            break;
                        }
                    }

                    scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, "Succes");
                }
                catch (Exception exception)
                {
                    string errorMessage = String.Format("Location '{0}', Soilprofile '{1}': {2}",
                                                        scenario.Location.Name, soilProfileProbability.SoilGeometryName, exception.Message);
                    PipingResults pipingResults = GetPipingResults(-1, pipingModelType, pipingCalculator);
                    scenario.SetPipingResults(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, pipingResults);
                    scenario.SetSafetyFactorPiping(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, -1);
                    scenario.Errors.Add("FAIL: " + errorMessage);
                    scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, "FAIL: " + errorMessage);
                }
            }
        }

        private PipingResults GetPipingResults(double? safetyFactor, PipingModelType pipingModelType, PipingCalculator pipingCalculator)
        {
            var pipingResults = new PipingResults();
            switch (pipingModelType)
            {
                case PipingModelType.Bligh:
                    pipingResults.BlighHCritical = pipingCalculator.HCritical;
                    pipingResults.BlighPipingFactor = safetyFactor;
                    break;
                case PipingModelType.Sellmeijer:
                    pipingResults.SellmeijerHCritical = pipingCalculator.HCritical;
                    pipingResults.SellmeijerPipingFactor = safetyFactor;
                    break;
                case PipingModelType.Sellmeijer2Forces:
                    pipingResults.Sellmeijer2ForcesHCritical = pipingCalculator.HCritical;
                    pipingResults.Sellmeijer2ForcesPipingFactor = safetyFactor;
                    break;
                case PipingModelType.Sellmeijer4Forces:
                    pipingResults.Sellmeijer4ForcesHCritical = pipingCalculator.HCritical;
                    pipingResults.Sellmeijer4ForcesPipingFactor = safetyFactor;
                    break;
                case PipingModelType.Wti2017:
                    pipingResults.Wti2017HCritical = pipingCalculator.HCritical;
                    pipingResults.Wti2017PipingFactor = safetyFactor;
                    break;
            }
            return pipingResults;
        }

        /// <summary>
        /// Ensures that the points on the surface line are never more than cDiff (0.5) apart.
        /// </summary>
        /// <param name="originalLine"></param>
        /// <returns></returns>
        private IEnumerable<GeometryPoint> GetCheckedSurfaceLine(IEnumerable<GeometryPoint> originalLine)
        {
            const double cDiff = 0.5;
            var newLine = new List<GeometryPoint>();
            double X = originalLine.First().X;
            foreach (var point in originalLine)
            {
                while (point.X > X + cDiff)
                {
                    var newPoint = new GeometryPoint(point)
                    {
                        X = X + cDiff
                    };
                    if (newPoint.X > newLine.Last().X)
                    {
                        newPoint.Z = newLine.Last().Z + ((newPoint.X - newLine.Last().X)/(point.X - newLine.Last().X))*
                                     (point.Z - newLine.Last().Z);
                        newLine.Add(newPoint);
                    }
                    X = newPoint.X;
                }
                newLine.Add(point);
            }
            return newLine;
        }

        private SurfaceLine2 DetermineNewSafeSurfaceLinePipingDeterministic(Scenario scenario, SoilGeometryProbability soilProfileProbability, PipingCalculator pipingCalculator, SurfaceLine2 surfaceLine, ref double? safetyFactor)
        {
            double orgShoulderLength = surfaceLine.DetermineShoulderWidth();
            double orgShoulderHeight = surfaceLine.DetermineShoulderHeight();
            double desiredShoulderLength = orgShoulderLength;
            double desiredShoulderHeight = orgShoulderHeight;
            double oldDesiredShoulderLength = orgShoulderLength;
            double oldDesiredShoulderHeight = orgShoulderHeight;

            GeometryPoint startSurfacePoint = surfaceLine.GetDikeToeInward();

            IEnumerable<GeometryPoint> relevantSurfacePointsList = from GeometryPoint point in surfaceLine.Geometry.Points
                                                                   where point.X >= startSurfacePoint.X
                                                                   orderby point.X ascending
                                                                   select point;
            relevantSurfacePointsList = GetCheckedSurfaceLine(relevantSurfacePointsList);
            int pointCount = 0;
            foreach (var point in relevantSurfacePointsList)
            {
                pointCount++;
                // Determine calculation filename to output piping calculation file
                //pipingCalculator.PipingCalculationDirectory = GetPipingCalculationBaseDirectory();
                //string fileNameCalculation =String.Format("Calc({0})_Loc({1})_Pro({2})_Pnt({3}))",
                //    pipingCalculator.CalculationModelIdentifier, scenario.Location.Name, soilProfileProbability.SoilProfile.Name, pointCount.ToString("d4")); ;
                //pipingCalculator.FilenameCalculation = Path.Combine(pipingCalculator.PipingCalculationDirectory, fileNameCalculation);

                // Calculate the piping design at the given point. This returns the required adaption (berm length and height) if any.
                var pipingDesign = pipingCalculator.CalculateDesignAtPoint(scenario.Location, surfaceLine,
                                                                           soilProfileProbability.SoilProfile,
                                                                           scenario.RiverLevel, point);
                if (pipingDesign != null)
                {
                    // Piping is an issue so adapt the surfaceline for it
                    desiredShoulderLength = pipingDesign.PipingLengthFromToe;
                    desiredShoulderLength = Math.Max(desiredShoulderLength, oldDesiredShoulderLength);
                    oldDesiredShoulderLength = desiredShoulderLength;
                    // shoulder height is height above surfacelevel!!
                    desiredShoulderHeight = pipingDesign.ShoulderHeightFromToe;
                    desiredShoulderHeight = Math.Max(desiredShoulderHeight, oldDesiredShoulderHeight);
                    oldDesiredShoulderHeight = desiredShoulderHeight;
                }
            }
            if (desiredShoulderLength > 0)
            {
                desiredShoulderLength = Math.Max(desiredShoulderLength, CMinimumShoulderWidth);
            }
            if (desiredShoulderLength > 0)
            {
                desiredShoulderHeight = Math.Max(desiredShoulderHeight, CMinimumShoulderElevation);
            }
            bool isNewShoulderSameAsOriginal = ((Math.Abs(desiredShoulderLength - orgShoulderLength) < CToleranceShoulderChanges) &&
                                                (Math.Abs(desiredShoulderHeight - orgShoulderHeight) < CToleranceShoulderChanges));
            if (isNewShoulderSameAsOriginal)
            {
                return surfaceLine;
            }
            else
            {
                // Adapt the surfaceline for the finally required shoulder dimensions.
                double maxShoulderLevel = CalculateMaximumShoulderLevel(surfaceLine, 1.0); // no limit to height of shoulder
                var surfaceLineShoulderAdapter = new SurfaceLineShoulderAdapter(surfaceLine, scenario.Location);
                surfaceLineShoulderAdapter.MaxShoulderLevel = maxShoulderLevel;
                SurfaceLine2 newSurfaceLine = surfaceLineShoulderAdapter.ConstructNewSurfaceLine(desiredShoulderLength, desiredShoulderHeight, true);

                safetyFactor = pipingCalculator.CalculatePipingFactor(scenario.Location, newSurfaceLine, soilProfileProbability.SoilProfile, scenario.RiverLevel);
                if (safetyFactor < scenario.RequiredSafetyFactorPiping)
                {
                    throw new DamFailureMechanismeCalculatorException("Deterministic Design: Piping is not safe yet.");
                }
                return newSurfaceLine;
            }
        }

        /// <summary>
        /// Create the initial geometry to be used to determine layers for embankment for design
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="stabilityCalculator"></param>
        /// <param name="soilProfileProbability"></param>
        /// <param name="soilGeometry2DName"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="initialgeometryFile"></param>
        private void CreateInitialGeometry(Scenario scenario, StabilityCalculator stabilityCalculator, SoilGeometryProbability soilProfileProbability, SurfaceLine2 surfaceLine, out String initialgeometryFile)
        {
            const int IterationIndex = -1;

            initialgeometryFile = StabilityCalculator.DetermineCalculationFilename(scenario.Location.Name, scenario.LocationScenarioID, soilProfileProbability.SoilGeometryName, IterationIndex);
            initialgeometryFile = initialgeometryFile + stabilityCalculator.GetFilenameExtension();
            initialgeometryFile = Path.Combine(stabilityCalculator.GetStabilityCalculationDirectory(), initialgeometryFile);
            double riverLevel = 0.5*(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).Z +
                                     surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).Z);
            string soilgeometry2DFilename = null;
            if (soilProfileProbability.SoilGeometry2DName != null)
            {
                soilgeometry2DFilename =
                    Path.GetFullPath(Path.Combine(DamProject.ProjectMap, Path.Combine(scenario.Location.MapForSoilGeometries2D, soilProfileProbability.SoilGeometry2DName)));
            }
            XDocument mstabXML = stabilityCalculator.CreateMStabXmlDoc(initialgeometryFile, scenario, soilProfileProbability.SoilProfile,
                                                                       soilgeometry2DFilename, riverLevel, null, surfaceLine);
            mstabXML.Save(initialgeometryFile + ".xml");
            var stabilityServiceAgent = new StabilityServiceAgent();
            stabilityServiceAgent.CreateProjectFile(mstabXML.ToString());
            if (!File.Exists(initialgeometryFile))
            {
                throw new DamFailureMechanismeCalculatorException("Initial geometry file (sti) is not created.");
            }
        }

        /// <summary>
        /// Redesigns the surface line for mechanism stability inside.
        /// </summary>
        /// <param name="scenario">The scenario.</param>
        /// <param name="requiredSafetyFactor">The safety factor.</param>
        private void RedesignSurfaceLineStabilityInside(Scenario scenario, double? requiredSafetyFactor)
        {
            Location location = scenario.Location;
            var modelParametersForPLLines = new ModelParametersForPLLines
            {
                PenetrationLength = location.PenetrationLength,
                DampingFactorPL3 = location.DampingFactorPL3,
                DampingFactorPL4 = location.DampingFactorPL4,
                PLLineCreationMethod = location.PLLineCreationMethod
            };
            using (var stabilityCalculator =
                new StabilityCalculator(damFailureMechanismeCalculationSpecification.FailureMechanismeParamatersMStab,
                                        programType,
                                        modelParametersForPLLines,
                                        location.TrafficLoad,
                                        location.MinimalCircleDepth,
                                        requiredSafetyFactor.Value,
                                        mstabProgramPath,
                                        slopeWProgramPath,
                                        location.GaugePLLines,
                                        location.Gauges,
                                        location.SoilbaseDB,
                                        location.SoilList,
                                        ProbabilisticType)
                {
                    SelectedStabilityKernelType = damFailureMechanismeCalculationSpecification.StabilityKernelType
                })
            {
                var validSoilProfileProbabilities = SelectStabilityProbabilities(scenario.Location.Segment.SoilProfileProbabilities);
                if (validSoilProfileProbabilities.Count == 0)
                {
                    throw new DamFailureMechanismeCalculatorException(String.Format("No stability profiles to calcutate for location {0}", scenario.Location.Name));
                }
                foreach (var soilProfileProbability in validSoilProfileProbabilities)
                {
                    // The calculation for UpliftVan will only will done if there is Uplift 
                    UpliftSituation? upliftSituation = scenario.GetStabilityUpliftSituation(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName);
                    if (IsCalculationRequired(upliftSituation))
                    {
                        switch (scenario.Location.StabilityDesignMethod)
                        {
                                // Bij WTIkernel is het soilprofile2D nog null in de soilProfileProbability en daarom gat Adapt Geometry fout.
// nagaan bij classic Delphi of dat daar ook leeg is. Daarnaast is het probleem groter. Je moet denk ik de geometrie geheel regenereren NADAT de surfaceline
// is aangepast. Dan moet de soilprofile2D goed gezet worden (en ook de sti + naam sti!?!?!?! uitzoeken!) op de nieuwe 2Dgeom. En die moet goed in de kernel terecht komen.
// Dus dat kan nog niet hier want de nieuwe surfaceline is nog niet bekend. Gebruik maken van de surf2dprof combiner!!
                            case StabilityDesignMethod.OptimizedSlopeAndShoulderAdaption:
                            {
                                RedesignCombinedSlopeAdaptionAndShoulderAdaption(scenario, stabilityCalculator, soilProfileProbability);
                                break;
                            }
                            case StabilityDesignMethod.SlopeAdaptionBeforeShoulderAdaption:
                            {
                                RedesignFirstSlopeAdaptionThenShoulderAdaption(scenario, stabilityCalculator, soilProfileProbability);
                                break;
                            }
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Redesign for Stability:
        /// Dependent on the exit point of the stability circle.
        /// If the exit point is in the slope, the slope will be adapted.
        /// If the exit point is not in the slope, a shoulder will be made, 
        /// or if the shoulder exists, the shoulder will be enlarged.
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="stabilityCalculator"></param>
        /// <param name="soilProfileProbability"></param>
        private void RedesignCombinedSlopeAdaptionAndShoulderAdaption(Scenario scenario, StabilityCalculator stabilityCalculator,
                                                                      SoilGeometryProbability soilProfileProbability)
        {
            Location location = scenario.Location;
            double requiredSafetyFactor = scenario.ModelFactors.RequiredSafetyFactorStabilityInnerSlope ??
                                          scenario.Location.ModelFactors.RequiredSafetyFactorStabilityInnerSlope.Value;
            double betaRequired = scenario.ModelFactors.RequiredProbabilityOfFailureStabilityInnerslope ??
                                  scenario.Location.ModelFactors.RequiredProbabilityOfFailureStabilityInnerslope.Value;
            const int maxRedesignIterations = 200;
            int iterationIndex = 1;
            string previousFilename;
            SurfaceLine2 orgSurfaceLine = scenario.Location.LocalXZSurfaceLine2;

            // Create the file with the initial geometry to be used to determine which layers have to be defined as dike embankment material 
            CreateInitialGeometry(scenario, stabilityCalculator, soilProfileProbability, orgSurfaceLine, out previousFilename);
            SurfaceLine2 surfaceLine = scenario.GetMostRecentSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName).FullDeepClone();

            var mstabDesignEmbankment = new MStabDesignEmbankment
            {
                EmbankmentMaterialname = location.DikeEmbankmentMaterial,
                PreviousGeometry2DFilename = previousFilename
            };

            try
            {
                double? beta;
                bool isRedesignRequired;
                double? exitPointXCoordinate;
                scenario.Location.AlignBoundaryPointsOfPL1LineWithAdaptedSurfaceLine(surfaceLine);
                stabilityCalculator.Calculate(scenario, soilProfileProbability.SoilProfile,
                                              GetFullSoilGeometry2DName(soilProfileProbability.SoilGeometry2DName),
                                              iterationIndex, mstabDesignEmbankment);
                mstabDesignEmbankment.PreviousGeometry2DFilename = stabilityCalculator.StabilityProjectFilename;
                mstabDesignEmbankment.EmbankmentMaterialname = location.ShoulderEmbankmentMaterial;
                MStabResults? mStabResults = scenario.GetMStabResults(soilProfileProbability.SoilProfile,
                                                                      soilProfileProbability.SoilGeometry2DName);
                double? safetyFactor = mStabResults.Value.zone1.safetyFactor;
                beta = scenario.GetFailureProbabilityStability(soilProfileProbability.SoilProfile,
                                                               soilProfileProbability.SoilGeometry2DName);
                isRedesignRequired = IsRedesignRequired(safetyFactor, requiredSafetyFactor, betaRequired, beta);
                exitPointXCoordinate = mStabResults.Value.zone1.circleSurfacePointRightXCoordinate;

                if (!isRedesignRequired && surfaceLine != null)
                {
                    // Set redesigned surfaceline to original, so in case no redesign is needed, the original surfaceline will be returned 
                    scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile,
                                                      soilProfileProbability.SoilGeometry2DName, surfaceLine);
                }
                else
                {
                    double maxFractionOfDikeHeightForShoulderHeight = scenario.Location.UseNewMaxHeightShoulderAsFraction ?
                                                                          scenario.Location.NewMaxHeightShoulderAsFraction : defaultMaxFractionOfDikeHeightForShoulderHeight;
                    double maxShoulderLevel = CalculateMaximumShoulderLevel(surfaceLine, maxFractionOfDikeHeightForShoulderHeight);
                    while (isRedesignRequired && surfaceLine != null)
                    {
                        iterationIndex++;
                        if (iterationIndex >= maxRedesignIterations)
                        {
                            throw new MaximumRedesignIterationsReachedException();
                        }

                        GeometryPoint limitPointForShoulderDesign = surfaceLine.GetLimitPointForShoulderDesign();
                        if (exitPointXCoordinate > limitPointForShoulderDesign.X)
                        {
                            // If exit point of circle is after the limitPointForShoulderDesign then enlarge the shoulder

                            // Determine new width and height for shoulder
                            double shoulderHeight;
                            double shoulderWidth;
                            DetermineNewShoulderWidthAndHeight(scenario.Location.StabilityShoulderGrowDeltaX,
                                                               scenario.Location.StabilityShoulderGrowSlope, surfaceLine, limitPointForShoulderDesign, out shoulderHeight, out shoulderWidth);

                            // Create new shoulder
                            var surfaceLineShoulderAdapter = new SurfaceLineShoulderAdapter(surfaceLine, scenario.Location);
                            surfaceLineShoulderAdapter.MaxShoulderLevel = maxShoulderLevel;
                            surfaceLine = surfaceLineShoulderAdapter.ConstructNewSurfaceLine(shoulderWidth, shoulderHeight, false);

                            var validationError = surfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                            if (validationError != null)
                            {
                                throw new SurfaceLineException(validationError.Text);
                            }

                            scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, surfaceLine);
                        }
                        else
                        {
                            // If exit point of circle is in the slope (inward) of the dike or the top of the shoulder then adapt slope
                            var surfaceLineSlopeAdapter = new SurfaceLineSlopeAdapter(surfaceLine, scenario.Location);
                            surfaceLine = surfaceLineSlopeAdapter.ConstructNewSurfaceLine(scenario.Location.StabilitySlopeAdaptionDeltaX);

                            var validationError = surfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                            if (validationError != null)
                            {
                                throw new SurfaceLineException(validationError.Text);
                            }
                            scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile,
                                                              soilProfileProbability.SoilGeometry2DName, surfaceLine);
                        }

                        scenario.Location.AlignBoundaryPointsOfPL1LineWithAdaptedSurfaceLine(surfaceLine);
                        stabilityCalculator.Calculate(scenario, soilProfileProbability.SoilProfile,
                                                      GetFullSoilGeometry2DName(soilProfileProbability.SoilGeometry2DName),
                                                      iterationIndex, mstabDesignEmbankment);

                        mStabResults = scenario.GetMStabResults(soilProfileProbability.SoilProfile,
                                                                soilProfileProbability.SoilGeometry2DName);
                        safetyFactor = mStabResults.Value.zone1.safetyFactor;
                        beta = scenario.GetFailureProbabilityStability(soilProfileProbability.SoilProfile,
                                                                       soilProfileProbability.SoilGeometry2DName);
                        isRedesignRequired = IsRedesignRequired(safetyFactor, requiredSafetyFactor, betaRequired, beta);
                        exitPointXCoordinate = mStabResults.Value.zone1.circleSurfacePointRightXCoordinate;
                    }
                }
                scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, "Succes");
            }
            catch (Exception exception)
            {
                string errorMessage = "FAIL: " + exception.Message;
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    errorMessage = errorMessage + ";" + innerException.Message;
                    innerException = innerException.InnerException;
                }
                scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, errorMessage);
                // Add message to log too for output tab.
                errorMessage = scenario.Location.Name + ", scenario " + scenario.LocationScenarioID + ", " + soilProfileProbability.SoilGeometryName + " " + errorMessage;
                stabilityCalculator.ErrorMessages.Add(new LogMessage(LogMessageType.Error, this, errorMessage));
                scenario.CalculationResult = CalculationResult.RunFailed;
                // Redesign not succesful, so no redesigned surfaceline will be returned 
                scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, null);
            }
            foreach (var errorMessage in stabilityCalculator.ErrorMessages)
            {
                errorMessages.Add(errorMessage);
            }
        }

        /// <summary>
        /// Determines the new height and width of the shoulder width
        /// </summary>
        /// <param name="shoulderGrowDeltaX">The shoulder grow delta X.</param>
        /// <param name="shoulderGrowSlope">The shoulder grow slope.</param>
        /// <param name="surfaceLine">The surface line.</param>
        /// <param name="limitPointForShoulderDesign">The limit point for shoulder design.</param>
        /// <param name="shoulderHeight">New height of the shoulder.</param>
        /// <param name="shoulderWidth">New width of the shoulder.</param>
        private static void DetermineNewShoulderWidthAndHeight(double shoulderGrowDeltaX, double shoulderGrowSlope,
                                                               SurfaceLine2 surfaceLine, GeometryPoint limitPointForShoulderDesign, out double shoulderHeight, out double shoulderWidth)
        {
            // Determine new shoulderpoint 
            var newShoulderPoint = new GeometryPoint()
            {
                X = limitPointForShoulderDesign.X + shoulderGrowDeltaX,
                Z = limitPointForShoulderDesign.Z + shoulderGrowDeltaX*shoulderGrowSlope
            };

            // Determine new shoulder width and height
            shoulderWidth = surfaceLine.DetermineShoulderLengthForGivenShoulderTopInside(newShoulderPoint);
            shoulderHeight = newShoulderPoint.Z - surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).Z;
        }

        /// <summary>
        /// Redesign for Stability:
        ///  - first do slope adaption (until required safety is reached)
        ///     * slope adaption is done in a range between startCoTangent and endCoTangent in steps of stepCoTangent
        ///     * and while the exit point of the slipcircle is in the slope (so currentCoTangent could be larger than endCoTangent
        ///  - than (if required safety is not reached) do shoulder adaption
        /// Note: Use same slope for shoulder as for dike.
        /// </summary>
        /// <param name="scenario"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="stabilityCalculator"></param>
        /// <param name="soilProfileProbability"></param>
        private void RedesignFirstSlopeAdaptionThenShoulderAdaption(Scenario scenario, StabilityCalculator stabilityCalculator,
                                                                    SoilGeometryProbability soilProfileProbability)
        {
            Location location = scenario.Location;
            double requiredSafetyFactor = scenario.ModelFactors.RequiredSafetyFactorStabilityInnerSlope ??
                                          scenario.Location.ModelFactors.RequiredSafetyFactorStabilityInnerSlope.Value;
            double betaRequired = scenario.ModelFactors.RequiredProbabilityOfFailureStabilityInnerslope ??
                                  scenario.Location.ModelFactors.RequiredProbabilityOfFailureStabilityInnerslope.Value;
            const int maxRedesignIterations = 200;
            int iterationIndex = 1;
            string previousFilename;
            SurfaceLine2 orgSurfaceLine = scenario.Location.LocalXZSurfaceLine2;

            // Create the file with the initial geometry to be used to determine which layers have to be defined as dike embankment material 
            CreateInitialGeometry(scenario, stabilityCalculator, soilProfileProbability, orgSurfaceLine, out previousFilename);
            SurfaceLine2 surfaceLine = scenario.GetMostRecentSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName).FullDeepClone();

            var mstabDesignEmbankment = new MStabDesignEmbankment
            {
                EmbankmentMaterialname = location.DikeEmbankmentMaterial,
                PreviousGeometry2DFilename = previousFilename
            };

            try
            {
                double? beta;
                bool isRedesignRequired;
                double? exitPointXCoordinate;
                scenario.Location.AlignBoundaryPointsOfPL1LineWithAdaptedSurfaceLine(surfaceLine);
                stabilityCalculator.Calculate(scenario, soilProfileProbability.SoilProfile,
                                              GetFullSoilGeometry2DName(soilProfileProbability.SoilGeometry2DName),
                                              iterationIndex, mstabDesignEmbankment);
                mstabDesignEmbankment.PreviousGeometry2DFilename = stabilityCalculator.StabilityProjectFilename;
                mstabDesignEmbankment.EmbankmentMaterialname = location.ShoulderEmbankmentMaterial;
                MStabResults? mStabResults = scenario.GetMStabResults(soilProfileProbability.SoilProfile,
                                                                      soilProfileProbability.SoilGeometry2DName);
                double? safetyFactor = mStabResults.Value.zone1.safetyFactor;
                beta = scenario.GetFailureProbabilityStability(soilProfileProbability.SoilProfile,
                                                               soilProfileProbability.SoilGeometry2DName);
                isRedesignRequired = IsRedesignRequired(safetyFactor, requiredSafetyFactor, betaRequired, beta);
                exitPointXCoordinate = mStabResults.Value.zone1.circleSurfacePointRightXCoordinate;

                if (!isRedesignRequired && surfaceLine != null)
                {
                    // Set redesigned surfaceline to original, so in case no redesign is needed, the original surfaceline will be returned 
                    scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile,
                                                      soilProfileProbability.SoilGeometry2DName, surfaceLine);
                }
                else
                {
                    double maxFractionOfDikeHeightForShoulderHeight = scenario.Location.UseNewMaxHeightShoulderAsFraction ?
                                                                          scenario.Location.NewMaxHeightShoulderAsFraction : defaultMaxFractionOfDikeHeightForShoulderHeight;
                    double maxShoulderLevel = CalculateMaximumShoulderLevel(surfaceLine, maxFractionOfDikeHeightForShoulderHeight);
                    while (isRedesignRequired && surfaceLine != null)
                    {
                        // First slope adaption
                        double startCoTangent = location.SlopeAdaptionStartCotangent;
                        double endCoTangent = location.SlopeAdaptionEndCotangent;
                        double stepCoTangent = location.SlopeAdaptionStepCotangent;
                        var orgCotangent = surfaceLine.GetCotangentOfInnerSlope();
                        double coTangent = startCoTangent;
                        double currentCoTangent = orgCotangent;

                        // Find out for which cotangent we want to start designing
                        while (coTangent <= orgCotangent)
                        {
                            coTangent += stepCoTangent;
                        }

                        // Design for slope adaption
                        GeometryPoint limitPointForShoulderDesign = surfaceLine.GetLimitPointForShoulderDesign();
                        while (isRedesignRequired && (coTangent < endCoTangent))
                        {
                            iterationIndex++;
                            if (iterationIndex >= maxRedesignIterations)
                            {
                                throw new MaximumRedesignIterationsReachedException();
                            }
                            var surfaceLineSlopeAdapter = new SurfaceLineSlopeAdapter(surfaceLine, scenario.Location);
                            // The parameter for ConstructNewSurfaceLineBySlope is the tangent of the slope, so use reciproce value
                            surfaceLine = surfaceLineSlopeAdapter.ConstructNewSurfaceLineBySlope(1/coTangent);
                            currentCoTangent = coTangent;

                            var validationError = surfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                            if (validationError != null)
                            {
                                throw new SurfaceLineException(validationError.Text);
                            }
                            scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile,
                                                              soilProfileProbability.SoilGeometry2DName, surfaceLine);

                            // Recalculate new surfaceline
                            scenario.Location.AlignBoundaryPointsOfPL1LineWithAdaptedSurfaceLine(surfaceLine);
                            stabilityCalculator.Calculate(scenario, soilProfileProbability.SoilProfile,
                                                          GetFullSoilGeometry2DName(soilProfileProbability.SoilGeometry2DName),
                                                          iterationIndex, mstabDesignEmbankment);
                            mStabResults = scenario.GetMStabResults(soilProfileProbability.SoilProfile,
                                                                    soilProfileProbability.SoilGeometry2DName);
                            safetyFactor = mStabResults.Value.zone1.safetyFactor;
                            exitPointXCoordinate = mStabResults.Value.zone1.circleSurfacePointRightXCoordinate;
                            beta = scenario.GetFailureProbabilityStability(soilProfileProbability.SoilProfile,
                                                                           soilProfileProbability.SoilGeometry2DName);
                            isRedesignRequired = IsRedesignRequired(safetyFactor, requiredSafetyFactor, betaRequired, beta);
                            limitPointForShoulderDesign = surfaceLine.GetLimitPointForShoulderDesign();

                            coTangent += stepCoTangent;
                        }

                        // Then shoulder adaption
                        while (isRedesignRequired && surfaceLine != null)
                        {
                            iterationIndex++;
                            if (iterationIndex >= maxRedesignIterations)
                            {
                                throw new MaximumRedesignIterationsReachedException();
                            }
                            
                            // Determine new width and height for shoulder
                            double shoulderHeight;
                            double shoulderWidth;
                            DetermineNewShoulderWidthAndHeight(scenario.Location.StabilityShoulderGrowDeltaX,
                                                               scenario.Location.StabilityShoulderGrowSlope, surfaceLine, limitPointForShoulderDesign, out shoulderHeight, out shoulderWidth);

                            // Create new shoulder
                            var surfaceLineShoulderAdapter = new SurfaceLineShoulderAdapter(surfaceLine, scenario.Location);
                            surfaceLineShoulderAdapter.MaxShoulderLevel = maxShoulderLevel;
                            surfaceLineShoulderAdapter.SlopeOfNewShoulder = currentCoTangent;
                            surfaceLine = surfaceLineShoulderAdapter.ConstructNewSurfaceLine(shoulderWidth, shoulderHeight, false);

                            var validationError = surfaceLine.Validate().FirstOrDefault(vr => vr.MessageType == ValidationResultType.Error);
                            if (validationError != null)
                            {
                                throw new SurfaceLineException(validationError.Text);
                            }
                            scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, surfaceLine);

                            // Recalculate new surfaceline
                            scenario.Location.AlignBoundaryPointsOfPL1LineWithAdaptedSurfaceLine(surfaceLine);
                            stabilityCalculator.Calculate(scenario, soilProfileProbability.SoilProfile,
                                                          GetFullSoilGeometry2DName(soilProfileProbability.SoilGeometry2DName),
                                                          iterationIndex, mstabDesignEmbankment);
                            mStabResults = scenario.GetMStabResults(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName);
                            safetyFactor = mStabResults.Value.zone1.safetyFactor;
                            beta = scenario.GetFailureProbabilityStability(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName);
                            isRedesignRequired = IsRedesignRequired(safetyFactor, requiredSafetyFactor, betaRequired, beta);
                            limitPointForShoulderDesign = surfaceLine.GetLimitPointForShoulderDesign();
                        }
                    }
                }
                scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, "Succes");
            }
            catch (Exception exception)
            {
                string errorMessage = "FAIL: " + exception.Message;
                Exception innerException = exception.InnerException;
                while (innerException != null)
                {
                    errorMessage = errorMessage + ";" + innerException.Message;
                    innerException = innerException.InnerException;
                }
                scenario.SetResultMessage(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, errorMessage);
                scenario.CalculationResult = CalculationResult.RunFailed;
                // Redesign not succesful, so no redesigned surfaceline will be returned 
                scenario.SetRedesignedSurfaceLine(soilProfileProbability.SoilProfile, soilProfileProbability.SoilGeometry2DName, null);
            }
            foreach (var errorMessage in stabilityCalculator.ErrorMessages)
            {
                errorMessages.Add(errorMessage);
            }
        }

        /// <summary>
        /// Redesigns the surface line flow slide.
        /// </summary>
        /// <param name="scenario">The scenario.</param>
        /// <param name="safetyFactor">The safety factor.</param>
        /// <param name="surfaceLine">The surface line.</param>
        private void RedesignSurfaceLineFlowSlide(Scenario scenario, double? safetyFactor, ref SurfaceLine2 surfaceLine)
        {
            throw new DamFailureMechanismeCalculatorException("Geometry adaption not supported for Flowslide");
        }

        /// <summary>
        /// Check if a stability redesign is required
        /// </summary>
        /// <param name="safetyFactor"></param>
        /// <param name="requiredSafetyFactor"></param>
        /// <param name="betaRequired"></param>
        /// <param name="beta"></param>
        /// <returns></returns>
        private bool IsRedesignRequired(double? safetyFactor, double requiredSafetyFactor, double betaRequired, double? beta)
        {
            bool isRedesignRequired;
            if (ProbabilisticType == ProbabilisticType.ProbabilisticFragility)
            {
                isRedesignRequired = beta > betaRequired;
            }
            else
            {
                isRedesignRequired = safetyFactor < requiredSafetyFactor;
            }
            return isRedesignRequired;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="pipingModelType"></param>
        /// <returns></returns>
        private PipingCalculator pipingCalculatorFactory(Location location, ModelParametersForPLLines modelParametersForPLLines, PipingModelType pipingModelType, double upliftCriterion, double requiredSafetyFactorPiping, PipingProbabilisticParameters? pipingProbabilisticParameters)
        {
            switch (pipingModelType)
            {
                case PipingModelType.Bligh:
                    return new PipingCalculatorBligh(modelParametersForPLLines,
                                                     requiredSafetyFactorPiping, location.GaugePLLines, location.Gauges, upliftCriterion);
                case PipingModelType.Sellmeijer:
                    PipingProbabilisticParameters? actualPipingProbabilisticParameters = null;
                    if ((ProbabilisticType == ProbabilisticType.Probabilistic) || (ProbabilisticType == ProbabilisticType.ProbabilisticFragility))
                    {
                        actualPipingProbabilisticParameters = pipingProbabilisticParameters;
                    }
                    return new PipingCalculatorSellmeijer(modelParametersForPLLines,
                                                          requiredSafetyFactorPiping, location.GaugePLLines, location.Gauges, actualPipingProbabilisticParameters, upliftCriterion);
                case PipingModelType.Sellmeijer2Forces:
                    return new PipingCalculatorSellmeijer2Forces(modelParametersForPLLines,
                                                                 requiredSafetyFactorPiping, location.GaugePLLines, location.Gauges, upliftCriterion);
                case PipingModelType.Sellmeijer4Forces:
                    return new PipingCalculatorSellmeijer4Forces(modelParametersForPLLines,
                                                                 requiredSafetyFactorPiping, location.GaugePLLines, location.Gauges, upliftCriterion);
                case PipingModelType.Wti2017:
                    return new PipingCalculatorWti2017(modelParametersForPLLines,
                                                       requiredSafetyFactorPiping, location.GaugePLLines, location.Gauges, upliftCriterion);
            }
            return null;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="soilgeometry2DName"></param>
        /// <returns></returns>
        private string GetFullSoilGeometry2DName(string soilGeometry2DName)
        {
            if (soilGeometry2DName == null)
            {
                return null;
            }
            else
            {
                string fullSoilGeometry2DName = Path.Combine(DamProject.ProjectMap, Path.Combine(MapForSoilGeometries2D, soilGeometry2DName));
                return fullSoilGeometry2DName;
            }
        }

        /// <summary>
        /// Get the directory where to create the piping project files
        /// </summary>
        /// <returns></returns>
        private static string GetPipingCalculationBaseDirectory()
        {
            return Path.Combine(Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location), "PipingCalculation\\");
        }

        /// <summary>
        /// Calculates the maximum level for the shoulder.
        /// </summary>
        /// <param name="surfaceLine">The surface line.</param>
        /// <param name="maxFractionOfDikeHeightForShoulderHeight">The fraction of dike height to determine maximimum shoulder height.</param>
        /// <returns></returns>
        private double CalculateMaximumShoulderLevel(SurfaceLine2 surfaceLine, double maxFractionOfDikeHeightForShoulderHeight)
        {
            var top = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).Z;
            var bottom = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).Z;
            if (top - bottom <= 0)
            {
                throw new SurfaceLineAdapterException(LocalizationManager.GetTranslatedText(this, "SurfaceLineShoulderAdapterMaxShoulderHeightError"));
            }
            double maxHeight = Math.Abs((top - bottom)*maxFractionOfDikeHeightForShoulderHeight);
            return bottom + maxHeight;
        }
    }
}