//-----------------------------------------------------------------------
// <copyright file="PipingCalculatorSellmeijer2Forces.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 2 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.Runtime.Serialization;
    

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

        [global::System.Serializable]
        public class PipingCalculatorSellmeijer2ForcesException : 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 PipingCalculatorSellmeijer2ForcesException() { }
            public PipingCalculatorSellmeijer2ForcesException(string message) : base(message) { }
            public PipingCalculatorSellmeijer2ForcesException(string message, Exception inner) : base(message, inner) { }
            protected PipingCalculatorSellmeijer2ForcesException(
              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 PipingCalculatorSellmeijer2ForcesException("An unexpected error occurred", e);
            }
        }

        static private double Factor(double whitesConstant, double beddingAngle)
        {
            double result = Math.PI / 3.0 * whitesConstant * CGammaPDividedByGammaW * Math.Tan(cDegreeToRadian * beddingAngle);
            return result;
        }

        static private double FRes(double whitesConstant, double beddingAngle)
        {
            const double cConst1 = 0.724789473684211;
            const double cConst2 = 0.50;
            double result = Factor(whitesConstant, beddingAngle) * Math.Pow((cConst2 / cConst1), 0.35);
            return result;
        }

        static private double FScale(double d70, double length, double permeability)
        {
            const double cConst3 = 207.921052631579;
            double cNuDividedByGravity = 0.00000133 / Physics.GravityConstant;
            double result = Math.Pow((d70 / cConst3), 0.39) * cConst3 * 0.000001 / Math.Pow(permeability * cNuDividedByGravity * length, 0.333);
            return result;
        }

        static private double FGeom(double dSandlayer, double length)
        {
            double result = 0.87 * Math.Pow(dSandlayer / length, 0.28 / (Math.Pow(dSandlayer / length, 2.8) - 1) + 0.04);
            return result;
        }

        static private double HcDividedByLNew(double whitesConstant, double beddingAngle, double d70, double length, double permeability, double dSandlayer)
        {
            double result = FRes(whitesConstant, beddingAngle) * FScale(d70, length, permeability) * FGeom(dSandlayer, length);
            return result;
        }

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

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

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

            HCritical = CalculateHCritical(
                pipingCalculatorInput.WhitesConstant,
                pipingCalculatorInput.BeddingAngle,
                Physics.FactorMeterToMicroMeter * aquiferLayerWithParameters.D70,
                pipingCalculatorInput.Length,
                aquiferLayerWithParameters.PermeabilityKx,
                aquiferLayerWithParameters.Height);
            // 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);

            // 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 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;
        }
    }
}