//-----------------------------------------------------------------------
// <copyright file="PipingCalculatorSellmeijer.cs" company="Deltares">
//  Copyright (c) 2009 Deltares. All rights reserved.
// </copyright>
// <author>R. Blankenburgh</author>
// <email>remko.blankenburgh@deltares.nl</email>
// <author>B. Faassen</author>
// <email>barry.faassen@deltares.nl</email>
// <date>07-07-2009</date>
// <summary>n.a.</summary>
//-----------------------------------------------------------------------

using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Probabilistic;

namespace Deltares.Dam.Data
{
    using System;
    using System.Linq;
    using System.Collections.Generic;
    using System.IO;
    using System.Diagnostics;
    using Deltares.Mathematics;
    using Deltares.Piping.Data;
    using Deltares.Standard;
    using System.Runtime.Serialization;
    using Deltares.Soilbase;

    public class PipingCalculatorSellmeijer : PipingCalculator
    {
        public PipingCalculatorSellmeijer(ModelParametersForPLLines modelParametersForPLLines, double requiredSafetyFactor,
            IList<GaugePLLine> gaugePLLines, IList<Gauge> gauges, PipingProbabilisticParameters? pipingProbabilisticParameters,
            double upliftCriterion) :
            base(modelParametersForPLLines, requiredSafetyFactor, gaugePLLines, gauges, pipingProbabilisticParameters, upliftCriterion)
        {
            CalculationModelIdentifier = "SellmeijerVNK";
        }

        [global::System.Serializable]
        public class PipingCalculatorSellmeijerException : ApplicationException
        {
            //
            // For guidelines regarding the creation of new exception types, see
            //    http://msdn.microsoft.com/library/default.asp?url=/library/en-us/cpgenref/html/cpconerrorraisinghandlingguidelines.asp
            // and
            //    http://msdn.microsoft.com/library/default.asp?url=/library/en-us/dncscol/html/csharp07192001.asp
            //

            public PipingCalculatorSellmeijerException() { }
            public PipingCalculatorSellmeijerException(string message) : base(message) { }
            public PipingCalculatorSellmeijerException(string message, Exception inner) : base(message, inner) { }
            protected PipingCalculatorSellmeijerException(
              SerializationInfo info,
              StreamingContext context)
                : base(info, context) { }
        }

        private double CalculatePipingFactorAtLevel(Location location, double waterLevel, double surfaceLevel)
        {
            var length = UpliftLocationAndResult.X - entryPoint.X;
            var dCoverLayer = DetermineHeightCoverLayer(surfaceLevel);

            var pipingModel2Calculation = CreateAndFillPipingModel2CalculationForDeterministic();
            pipingModel2Calculation.IsHeadDropCalculation = true;// Constant
            pipingModel2Calculation.SeepageLength = length; // L Calculated
            pipingModel2Calculation.CrackLength = dCoverLayer; // Distance upper sandlayer to surfaceline, "Do" in DGSWebPiping
            // Cracklength = FluidisationLength = VerticalSeepageLength

            // Reference level is highest value of surfaceLevel or PolderLevel
            // Uit TR Zandmeevoerende wellen (1999): "Het verval dH is gelijk aan het verschil tussen buitenwaterstand (het ontwerppeil(OP)) 
            // bij zeedijken en de maatgevende hoogwaterstand (MHW bij rivierdijken) en de waterstand binnendijks ter plaatse van het uittredepunt, 
            // rekening houdend met zeespiegelrijzing etc.(zie paragraaf 3.7.2). In dien ter plaatse van het uittreepunt of de opbarstlocatie 
            // geen vrije waterstand heerst kan gerekend worden met het maaiveldniveau, rekening houdend met eventuele maaiveld daling (zie paragraaf 3.7.2)."
            var referenceLevel = Math.Max(location.PolderLevel, surfaceLevel);

            // Ouput piping parameters and dam parameters to project file
            var fileName = "";
            if (FilenameCalculation != null)
            {
                if (FilenameCalculation != "")
                {
                    string directoryName = Path.GetDirectoryName(FilenameCalculation);
                    if (!String.IsNullOrEmpty(directoryName) && !Directory.Exists(directoryName))
                    {
                        Directory.CreateDirectory(directoryName);

                    }
                    fileName = String.Format("{0}{1}", FilenameCalculation, PipingFilenameExtension);
                    Console.WriteLine(fileName);
                    File.WriteAllLines(fileName, pipingModel2Calculation.ParametersDeterministicToStrings().ToArray());

                    // Write Dam parameters to calculation file
                    var damParameterStrings = new List<string>
                                                  {
                                                      DamParametersSectionId,
                                                      String.Format("PolderLevel={0}", location.PolderLevel),
                                                      String.Format("SurfaceTopLevel={0}", UpliftLocationAndResult.Z),
                                                      String.Format("ReferenceLevel={0}", referenceLevel),
                                                      String.Format("WaterLevel={0}", waterLevel),
                                                      ""
                                                  };
                    File.AppendAllText(fileName, string.Join(Environment.NewLine, damParameterStrings.ToArray()));
                }
            }
            pipingModel2Calculation.CalculateHeadDropPC2();

            HCritical = pipingModel2Calculation.HeadDrop;
            // The Leidraad says that 0.3 * d can be substracted from hActual
            // In the piping Sellmeijer model (PipingModel2 for VNK) 0.3 * d  is added to the hCritical
            // So we do not substract 0.3 * d from hActual (written by Tom The and Erik Vastenburg)
            var hActual = waterLevel - referenceLevel;
            var pipingFactor = cDefaultMaxReturnValue;
            // If actual headdrop is almost zero or smaller than zero then piping will be no problem 
            // The pipingfactor will be set to cDefaultMaxReturnValue
            if (hActual > cToleranceHead)
            {
                pipingFactor = HCritical / hActual;
            }

            // Write results to calculation file
            if (fileName != "")
            {
                var resultStrings = new List<string>
                                        {
                                            ResultsSectionId,
                                            String.Format("HCritical={0}", HCritical),
                                            String.Format("HActual={0}", hActual),
                                            String.Format("PipingFactor={0}", pipingFactor),
                                            ""
                                        };
                File.AppendAllText(fileName, string.Join(Environment.NewLine, resultStrings.ToArray()));
            }
            return pipingFactor;            
        }

        public override PipingDesign CalculateDesignAtPoint(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, double waterLevel, GeometryPoint point)
        {
            // call base class for initializations, this will also calculate the uplift
            var pipingDesign = base.CalculateDesignAtPoint(location, surfaceLine, soilProfile, waterLevel, point);
            // if there is nu uplift, then there is no piping so return null
            if ((UpliftLocationAndResult != null) && (UpliftLocationAndResult.UpliftFactor != null) && (UpliftLocationAndResult.UpliftFactor.Value < upliftCriterion))
            {
                // Calculate the piping safety factor at the given point
                var sp = CalculatePipingFactorAtLevel(location, waterLevel, point.Z);

                // If too low, then determine required height and length (from uplift)
                if (sp < requiredSafetyFactor)
                {
                    // Finally, determine the required shoulderheight
                    double currentShoulderHeight = UpliftLocationAndResult.Z -
                                                   surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).Z;
                    pipingDesign.ShoulderHeightFromToe = currentShoulderHeight + CalculateExtraShoulderHeight();
                    pipingDesign.PipingLengthFromToe = UpliftLocationAndResult.X - surfaceLine.GetDikeToeInward().X;
                    return pipingDesign;
                }
            }
            return null;
        }

        /// <summary>
        /// Calculate the reliability index
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="headPL2"></param>
        /// <param name="polderLevel"></param>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        public override double? CalculateReliabilityIndex(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile,
            ProbabilisticStruct waterLevel)
        {
            base.CalculateReliabilityIndex(location, surfaceLine, soilProfile, waterLevel);
            ThrowIfNoPLLinesDefined();
            ThrowHelper.ThrowIfArgumentNull(surfaceLine, StringResourceNames.SurfaceLineNotAssigned);
            ThrowHelper.ThrowIfArgumentNull(soilProfile, StringResourceNames.SoilProfileNotAssigned);

            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;
            double? reliabilityIndex = null;
            double probabilityOfFailure = 0.0;
            int pointIndex = 0;
            foreach (GeometryPoint surfacePoint in relevantSurfacePointsList)
            {
                pointIndex++;
                UpliftLocationDeterminator upliftLocationDeterminator = new UpliftLocationDeterminator
                {
                    PLLines = this.PLLines,
                    SoilProfile = soilProfile,
                    SurfaceLine = surfaceLine,
                    XSoilGeometry2DOrigin = location.XSoilGeometry2DOrigin
                };
                UpliftLocationAndResult = upliftLocationDeterminator.GetUpliftFactorAtPoint(surfacePoint);
                double reliabilityIndexUplift = UpliftFactorToBetaCalculator.ComputeBetaFromUpliftFactor(UpliftLocationAndResult.UpliftFactor.Value);

                try
                {
                    double seepageLength = UpliftLocationAndResult.X - entryPoint.X;
                    double topLevel = SoilProfile.GetLayerWithId(UpliftLocationAndResult.LayerWhereUpliftOccuresId).TopLevel;
                    double d = UpliftLocationAndResult.Z - topLevel;

                    var pipingModel2Calculation = CreateAndFillPipingModel2CalculationForProbabilistic();

                    ProbabilisticStruct probWaterlevel = waterLevel;
                    // waterLevel is specified in NAP. Do not forget to substract the reference level, because we have to calculate with delta H !!
                    double referenceLevel = Math.Max(location.PolderLevel, topLevel);
                    probWaterlevel.Mean = Math.Max(0.3, waterLevel.Mean - referenceLevel);
                    pipingModel2Calculation.ProbExternalWaterlevel = probWaterlevel;

                    ProbabilisticStruct probSeepageLength = new ProbabilisticStruct(seepageLength, 0.0, (int)DistributionType.Deterministic, false);
                    pipingModel2Calculation.ProbSeepageLength = probSeepageLength;

                    pipingModel2Calculation.IsHeadDropCalculation = true;// Constant
                    pipingModel2Calculation.SeepageLength = seepageLength; // L Calculated
                    pipingModel2Calculation.CrackLength = d; // Distance upper sandlayer to surfaceline, "Do" in DGSWebPiping
                    // Cracklength = FluidisationLength = VerticalSeepageLength

                    // Write project file
                    if (this.FilenameCalculation != "")
                    {
                        File.WriteAllLines(String.Format("{0}_Pnt{1}{2}", this.FilenameCalculation, pointIndex, PipingFilenameExtension), pipingModel2Calculation.ParametersProbabilisticToStrings().ToArray());
                    }
                    double reliabilityIndexPiping = pipingModel2Calculation.CalculateReliabilityIndex();
                    probabilityOfFailure = Math.Max(probabilityOfFailure,
                                                    Probabilistic.Probabilistic.NormalDistribution(-reliabilityIndexUplift) *
                                                    Probabilistic.Probabilistic.NormalDistribution(-reliabilityIndexPiping));
                }
                catch (Exception e)
                {
                    throw new DamFailureMechanismeCalculatorException("Couldn't calculate the reliability using the Sellmeijer calculator", e);
                }
            }
            reliabilityIndex = -Probabilistic.Probabilistic.NormalDistrInverse(probabilityOfFailure);

            return reliabilityIndex;
        }

        /// <summary>
        /// Calculate the probability of failure
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="headPL2"></param>
        /// <param name="polderLevel"></param>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        public override double? CalculatePipingFailureProbability(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile,
            ProbabilisticStruct waterLevel)
        {
            double? reliabilityIndex;
            reliabilityIndex = CalculateReliabilityIndex(location, surfaceLine, soilProfile, waterLevel);
            if (reliabilityIndex != null)
            {
                return Probabilistic.Probabilistic.NormalDistribution(-reliabilityIndex.Value);
            }
            else
            {
                return null;
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        public override double? CalculatePipingFactor(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, double waterLevel)
        {
            base.CalculatePipingFactor(location, surfaceLine, soilProfile, waterLevel);
            try
            {
                if (UpliftLocationAndResult != null)
                {
                    return CalculatePipingFactorAtLevel(location, waterLevel,
                                                        surfaceLine.Geometry.GetZAtX(UpliftLocationAndResult.X));                    
                }
                return cDefaultMaxReturnValue;
            }
            catch (Exception e)
            {
                throw new DamFailureMechanismeCalculatorException("Couldn't calculate the head drop using the Sellmeijer calculator", e);
            }
        }

        /// <summary>
        /// Create PipingModel2Calculation object with deterministic data filled from model
        /// </summary>
        /// <returns></returns>
        private PipingModel2Calculation CreateAndFillPipingModel2CalculationForDeterministic()
        {
            var pipingModel2Calculation = new PipingModel2Calculation();

            PipingCalculatorInput pipingCalculatorInput = DefinePipingInputFromModel();

            pipingModel2Calculation.HeadDrop = 0.0;// Not used, to be calculated
            pipingModel2Calculation.PipingCommonData.BeddingAngle = pipingCalculatorInput.BeddingAngle; // From soilmaterials of upper sandlayer
            pipingModel2Calculation.PipingCommonData.FluidisationGradient = 0.3; //Constant
            pipingModel2Calculation.PipingCommonData.WaterUnitWeight = Physics.UnitWeightOfwater; //Constant
            pipingModel2Calculation.PipingCommonData.WaterViscosity = pipingCalculatorInput.WaterViscosity; //Constant
            pipingModel2Calculation.PipingCommonData.WhitesConstant = pipingCalculatorInput.WhitesConstant; // From soilmaterials of upper sandlayer
            pipingModel2Calculation.PipingCommonData.SafetyFactor = 1.20; //Constant toetsen/ontwerpen = 1.2 en calamiteit = 1.0
            pipingModel2Calculation.PipingCommonData.IsAdjustHeadDrop = true; // Constant

            // Soil parameters
            pipingModel2Calculation.PipingCommonData.ParticleUnitWeight = 26.50; // Constant
            pipingModel2Calculation.Permeability1 = pipingCalculatorInput.PermeabilityInBetweenAquiferlayer; // Kx From upper sandlayer
            pipingModel2Calculation.Permeability2 = pipingCalculatorInput.PermeabilityBottomAquiferlayer; // Kx From bottom sandlayer
            pipingModel2Calculation.Permeability3 = pipingCalculatorInput.PermeabilityInBetweenAquiferlayer; // Kx From upper sandlayer
            pipingModel2Calculation.ParticleDiameter = pipingCalculatorInput.D70;  // D70 D60 From soilmaterials of upper sandlayer

            // Geometry parameters
            pipingModel2Calculation.Height1 = pipingCalculatorInput.DInBetweenAquiferlayer; // Thickness upper sandlayer 
            pipingModel2Calculation.Height2 = pipingCalculatorInput.DBottomAquiferlayer; // Thickness bottom sandlayer 
            return pipingModel2Calculation;
            //
            // particleunitweight = gamma_p can be calculated as follows
            // (gamma_sat - gamma_water) = (gamma_p - gamma_water) / (1 - n); for now 26.50 is good enough
            //
            // limit heights: see TPC2Preprocessing.AdjustParameters()
            // CD1MinSeepageLengthFactor = 0.02; = D1 / Length 
            // CD2MinSeepageLengthFactor = 0.03; = D2 / Length
            // CD1MaxSeepageLengthFactor = 1.4; 
            // CD2MaxSeepageLengthFactor = 1.4; 
            //
            // Kx or Ky: Sel: Use Kx
        }

        /// <summary>
        /// Create PipingModel2Calculation object with probabilistic data filled from model
        /// </summary>
        /// <returns></returns>
        private PipingModel2Calculation CreateAndFillPipingModel2CalculationForProbabilistic()
        {
            var pipingModel2Calculation = new PipingModel2Calculation();
            PipingCalculatorInput pipingCalculatorInput = DefinePipingInputFromModel();
            pipingModel2Calculation.HeadDrop = 0.0;// Not used
            pipingModel2Calculation.PipingCommonData.FluidisationGradient = CDefaultFluidisationGradient; //Constant
            pipingModel2Calculation.PipingCommonData.WaterUnitWeight = Physics.UnitWeightOfwater; //Constant
            pipingModel2Calculation.PipingCommonData.WaterViscosity = pipingCalculatorInput.WaterViscosity; //Constant
            pipingModel2Calculation.PipingCommonData.SafetyFactor = 1.20; //Constant toetsen/ontwerpen = 1.2 en calamiteit = 1.0
            pipingModel2Calculation.PipingCommonData.IsAdjustHeadDrop = true;
            pipingModel2Calculation.CrackLength = 0.1;
            pipingModel2Calculation.IsHeadDropCalculation = true;

            // Soil parameters
            pipingModel2Calculation.PipingCommonData.ProbBeddingAngle = pipingCalculatorInput.ProbBeddingAngle;
            pipingModel2Calculation.PipingCommonData.ProbParticleD70 = pipingCalculatorInput.ProbD70;
            pipingModel2Calculation.PipingCommonData.ProbParticleUnitWeight = new ProbabilisticStruct(26.50, 0.0, (int)DistributionType.Deterministic, false);
            pipingModel2Calculation.PipingCommonData.ProbWhitesConstant = pipingCalculatorInput.ProbWhitesConstant;

            // Geometry parameters
            if (this.PipingProbabilisticParameters != null)
            {
                ProbabilisticStruct probabilisticParameter;
                probabilisticParameter = new ProbabilisticStruct(pipingCalculatorInput.DInBetweenAquiferlayer, this.PipingProbabilisticParameters.Value.LayerHeightDeviation, (int)this.PipingProbabilisticParameters.Value.LayerHeightDistribution, true);
                pipingModel2Calculation.ProbLayer1Thickness = probabilisticParameter;
                probabilisticParameter = new ProbabilisticStruct(pipingCalculatorInput.DBottomAquiferlayer, this.PipingProbabilisticParameters.Value.LayerHeightDeviation, (int)this.PipingProbabilisticParameters.Value.LayerHeightDistribution, true);
                pipingModel2Calculation.ProbLayer2Thickness = probabilisticParameter;
            }
            else
            {
                ProbabilisticStruct probabilisticParameter;
                probabilisticParameter = new ProbabilisticStruct(pipingCalculatorInput.DInBetweenAquiferlayer, 0, 0, false);
                pipingModel2Calculation.ProbLayer1Thickness = probabilisticParameter;
                probabilisticParameter = new ProbabilisticStruct(pipingCalculatorInput.DBottomAquiferlayer, 0, 0, false);
                pipingModel2Calculation.ProbLayer2Thickness = probabilisticParameter;
            }
            pipingModel2Calculation.ProbSoil1Permeability = pipingCalculatorInput.ProbPermeabilityInBetweenAquiferlayer;
            pipingModel2Calculation.ProbSoil2Permeability = pipingCalculatorInput.ProbPermeabilityBottomAquiferlayer;
            pipingModel2Calculation.ProbSoil3Permeability = pipingCalculatorInput.ProbPermeabilityInBetweenAquiferlayer;

            return pipingModel2Calculation;
        }
    }
}