//-----------------------------------------------------------------------
// <copyright file="PipingCalculatorSellmeijer4Forces.cs" company="Deltares">
//  Copyright (c) 2010 Deltares. All rights reserved.
// </copyright>
// <author>B.S.T.I.M. The</author>
// <email>tom.the@deltares.nl</email>
// <date>13-7-2010</date>
// <summary>Calculates the piping factor according to Sellmeijer 4 forces.</summary>
//-----------------------------------------------------------------------

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

namespace Deltares.Dam.Data
{
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Runtime.Serialization;
    using Deltares.Soilbase;

    public class PipingCalculatorSellmeijer4Forces : PipingCalculator
    {
        public PipingCalculatorSellmeijer4Forces(ModelParametersForPLLines modelParametersForPLLines, double requiredSafetyFactor,
            IList<GaugePLLine> gaugePLLines, IList<Gauge> gauges, double requiredUpliftFactor) :
            base(modelParametersForPLLines, requiredSafetyFactor, gaugePLLines, gauges, null,
            requiredUpliftFactor)
        {
        }

        private const double cBligh = 18.0;

        [global::System.Serializable]
        public class PipingCalculatorSellmeijer4ForcesException : 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 PipingCalculatorSellmeijer4ForcesException() { }
            public PipingCalculatorSellmeijer4ForcesException(string message) : base(message) { }
            public PipingCalculatorSellmeijer4ForcesException(string message, Exception inner) : base(message, inner) { }
            protected PipingCalculatorSellmeijer4ForcesException(
              SerializationInfo info,
              StreamingContext context)
                : base(info, context) { }
        }

        override public 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 PipingCalculatorSellmeijer4ForcesException("An unexpected error occurred", e);
            }
        }

        static public double TermLnC(double whitesConstant, double beddingAngle, double d70, double length, double permeability, double waterViscosity)
        {
            double result = 0.68 - 0.1 * Math.Log(TermC(whitesConstant, beddingAngle, d70, length, permeability, waterViscosity), Math.E);
            return result;
        }
        static public double TermC(double whitesConstant, double beddingAngle, double d70, double length, double permeability, double waterViscosity)
        {
            double kIntrinsic = KIntrinisc(permeability, waterViscosity);
            double term = Math.Pow((1 / (kIntrinsic * length)), (1.0 / 3.0));
            double result = whitesConstant * (d70 * 1.0e-6) * term;
            return result;
        }

        static public double TermAlpha(double length, double dSandlayer)
        {
            double result = Math.Pow((dSandlayer / length), 0.28 / ((Math.Pow(dSandlayer / length, 2.8)) - 1));
            return result;
        }

        static private double HcDividedByLNew(double whitesConstant, double beddingAngle, double d70, double length, double permeability, double dSandlayer, double waterViscosity)
        {
            return TermAlpha(length, dSandlayer) * TermC(whitesConstant, beddingAngle, d70, length, permeability, waterViscosity) * 
                CGammaPDividedByGammaW * Math.Tan(cDegreeToRadian * beddingAngle) * 
                TermLnC(whitesConstant, beddingAngle, d70, length, permeability, waterViscosity);
             
        }

        static public double CalculateHCritical(double whitesConstant, double beddingAngle, double d70, double length, double permeability, double dSandlayer, double waterViscosity)
        {
            double result = HcDividedByLNew(whitesConstant, beddingAngle, d70, length, permeability, dSandlayer,
                                             waterViscosity) * length;
            return result;
        }

        private double CalculatePipingFactorAtLevel(Location location, double waterLevel, double surfaceLevel)
        {
            PipingCalculatorInput pipingCalculatorInput = DefinePipingInputFromModel();

            var dCoverLayer = DetermineHeightCoverLayer(surfaceLevel);
            var soilProfile1DAquiferLayerCombiner = new SoilProfile1DAquiferLayerCombiner(SoilProfile);
            AquiferLayerWithParameters aquiferLayerWithParameters = soilProfile1DAquiferLayerCombiner.CombineLayers(UpliftLocationAndResult.LayerWhereUpliftOccuresId);

            this.HCritical = CalculateHCritical(
                pipingCalculatorInput.WhitesConstant,
                pipingCalculatorInput.BeddingAngle,
                Physics.FactorMeterToMicroMeter * aquiferLayerWithParameters.D70,
                pipingCalculatorInput.Length,
                aquiferLayerWithParameters.PermeabilityKx,
                aquiferLayerWithParameters.Height,
                pipingCalculatorInput.WaterViscosity);

            // 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)."
            double referenceLevel = Math.Max(location.PolderLevel, surfaceLevel);

            // After consulting Erik Vastenburg, Jan Blinde, Andre Koelewijn, Alexander van Duinen en Huub de Bruijn
            // it was decided that the fluidisation correction shopuld be applied on the hActual and not the hCritical.
            double hActual = waterLevel - referenceLevel - (CDefaultFluidisationGradient * dCoverLayer);
            double pipingFactor = DeterminePipingFactorWithUpperLimit(hActual, HCritical);
            return pipingFactor;            
        }



        /// <summary>
        /// Determines the required shoulder (height and length as seen from the original dike toe) 
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="waterLevel"></param>
        /// <param name="point"></param>
        /// <returns>The required dimensions when needed else null</returns>
        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)
            {
                // Calculate the piping safety factor using the level of 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;
        }

    }
}