///
/// Copyright (c) 2013-2018 Sensus Slovensko a.s.
///
using System;
using System.Collections.Generic;
using log4net;

namespace Config
{
    public static class Formulas
    {
		private static readonly ILog log = LogManager.GetLogger(typeof(Formulas));

        ///
        /// Tables to calculate specific enthalpy
        ///
        private static readonly int[] Ii;
        private static readonly int[] Ji;
        private static readonly double[] ni;

		private static readonly double[] Di;

        /// <summary> 
		/// Constructor
		/// </summary>
        static Formulas()
        {
            ///
            /// Initialize tables to calculate specific enthalpies
            ///
            Ii = new int[34] { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32 };
            Ji = new int[34] { -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41 };
            ni = new double[34] {
                 0.14632971213167,      /// 1
                -0.84548187169114,      /// 2
                -0.37563603672040E1,    /// 3
                 0.33855169168385E1,    /// 4
                -0.95791963387872,      /// 5
                 0.15772038513228,      /// 6
                -0.16616417199501E-1,   /// 7
                 0.81214629983568E-3,   /// 8
                 0.28319080123804E-3,   /// 9
                -0.60706301565874E-3,   /// 10
                -0.18990068218419E-1,   /// 11
                -0.32529748770505E-1,   /// 12
                -0.21841717175414E-1,   /// 13
                -0.52838357969930E-4,   /// 14
                -0.47184321073267E-3,   /// 15
                -0.30001780793026E-3,   /// 16
                 0.47661393906987E-4,   /// 17
                -0.44141845330846E-5,   /// 18
                -0.72694996297594E-15,  /// 19
                -0.31679644845054E-4,   /// 20
                -0.28270797985312E-5,   /// 21
                -0.85205128120103E-9,   /// 22
                -0.22425281908000E-5,   /// 23
                -0.65171222895601E-6,   /// 24
                -0.14341729937924E-12,  /// 25
                -0.40516996860117E-6,   /// 26
                -0.12734301741641E-8,   /// 27
                -0.17424871230634E-9,   /// 28
                -0.68762131295531E-18,  /// 29
                 0.14478307828521E-19,  /// 30
                 0.26335781662795E-22,  /// 31
                -0.11947622640071E-22,  /// 32
                 0.18228094581404E-23,  /// 33
                -0.93537087292458E-25,  /// 34
            };

			///
			/// Initialize a table to calculate temperature of a platinum thermometer from resistance
			///
			Di = new double[]
			{
				 439.932854,
				 472.418020,
				 37.684494,
				 7.472018,
				 2.920828,
				 0.005184,
				-0.963864,
				-0.188732,
				 0.191203,
				 0.049025,
			};
        }

        /// <summary>
        /// Calculate density of distilled water from temperature
        /// </summary>
        /// <param name="t">ITS-90 temperature in [°C]</param>
        /// <returns>Density in [kg/m3]</returns>
        public static double DistilledWaterDensityFromTemp(double t)
        {
            if (t <= 40)
            {
                const double c0 = 999.839564;
                const double c1 = 0.067998613;
                const double c2 = -0.0091101468;
                const double c3 = 0.00010058299;
                const double c4 = -0.0000011275659;
                const double c5 = 6.5985371e-09;

                return ((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0;
            }
            else
            {
                const double a0 = 9.9983952E2;
                const double a1 = 1.6952577E1;
                const double a2 = -7.9905127E-3;
                const double a3 = -4.6241757E-5;
                const double a4 = 1.0584601E-7;
                const double a5 = -2.8103006E-10;
                const double b = 1.6887236E-2;

                return (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t + a0) / (1.0 + b * t);
            }
        }

		/// <summary>
		/// Calculate density of distilled water from temperature (obsolete)
		/// </summary>
		/// <param name="t">IPTS-68 temperature in [°C]</param>
		/// <returns>Density in [kg/m3]</returns>
        public static double DistilledWaterDensityFromTempIPTS68(double t)
        {
            const double a0 = 999.842594;
            const double a1 = 0.06793952;
            const double a2 = -0.009095290;
            const double a3 = 0.0001001685;
            const double a4 = -0.000001120083;
            const double a5 = 6.536332e-09;

            return ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t + a0;
        }

		/// <summary>
		/// Calculate density by comparing calculated data and data from a certificate
		/// </summary>
		/// <param name="realDensity">Density from a certificate in [kg/m3]</param>
		/// <param name="atTemperature">Temperature from a certificate in [°C]</param>
		/// <returns>Density correction in [kg/m3]</returns>
		public static double DensityCorrection(double realDensity, double atTemperature)
		{
			/// Calculated data
			double calculatedDensity = DistilledWaterDensityFromTemp(atTemperature);

			return realDensity - calculatedDensity;
		}

		/// <summary>
		/// Calculate corrected (real) water density from temperature
		/// </summary>
		/// <param name="t">Temperature in [°C]</param>
		/// <returns>Density in [kg/m3]</returns>
		public static double WaterDensityFromTemp(double t, double realDensity, double atTemperature)
		{
            return DistilledWaterDensityFromTemp(t) + DensityCorrection(realDensity, atTemperature);
		}

        /// <summary>
        /// Calculate corrected (real) water density from temperature
        /// </summary>
        /// <param name="t">Temperature in [°C]</param>
        /// <returns>Density in [kg/m3]</returns>
        public static double WaterDensityFromTempPress(double temp, double pressure, double realDensity, double atTemperature)
        {
            double x0 = 5.08821E-10;
            double x1 = 1.2639418;
            double x2 = 0.2660269;
            double x3 = 0.3734838;
            double x4 = 2.0205242;
            double theta = temp / 100.0;
            double B = x0 * (((x3 * theta + x2) * theta + x1) * theta + 1) / (1 + x4 * theta);

            return WaterDensityFromTemp(temp, realDensity, atTemperature) * (1 + B * Config.Units.ConvertTo(Config.Unit.Pa, pressure));
        }


		public static float AirDensityFromAmbientVales(float tempC, float pressureBar, float humiPct)
		{
			double pressurePa = 100000.0 * (double)pressureBar;		/// [Pa]
			double tempKelvin = 273.15 + (double)tempC;
			double coef1 = 1.2811805 / 10000.0 * tempKelvin * tempKelvin
						 - 1.950987 / 100.0 * tempKelvin
						 + 34.04926034
						 - 6.353631 * 1000.0 / tempKelvin;
			double coef3 = humiPct / 100.0 * System.Math.Exp(coef1) / pressurePa;
			double airDensityKgm3 = 0.00348353 * pressurePa * (1.0 - 0.378 * coef3) / tempKelvin;	/// kg/m3
			return (float)airDensityKgm3;
		}

		/// <summary>
		/// Convert 'pulses' to 'volume', prevent division by zero
		/// </summary>
		public static double VolumeFromPulses(int pulses, double pulsesPerLiter)
		{
			if (pulsesPerLiter <= double.Epsilon) return 0;
			return Convert.ToDouble(pulses) / pulsesPerLiter;
		}

		/// <summary>
		/// Calculate the error in % from 'measured' and 'true' volume, prevent division by zero
		/// </summary>
		public static double ErrorFromVolumes(double measuredVolume, double trueVolume)
		{
			if (trueVolume <= float.Epsilon)
			{
				if (measuredVolume <= float.Epsilon)
				{
					log.WarnFormat("ErrorFromVolumes({0},{1}) returns {2}", measuredVolume, trueVolume, -100.0);
					return -100.0;
				}

				log.WarnFormat("ErrorFromVolumes({0},{1}) returns {2}", measuredVolume, trueVolume, 99.0);
				return 99.0;
			}

			double error = 100.0 * (measuredVolume - trueVolume) / trueVolume;
			log.InfoFormat("ErrorFromVolumes({0},{1}) returns {2}", measuredVolume, trueVolume, error);
			return error;
		}

		/// <summary>
		/// Calculates corrected value form a list of corrections by interpolation.
		/// It is assumed that values in the list 'corrections' are sorted.
		/// </summary>
		/// <param name="rawMeasurement">Raw uncorrected value</param>
		/// <param name="corrections">Sorted (value, correction) pairs</param>
		/// <returns>Corrected value</returns>
		public static double CorrectedValue(double rawValue, IList<Config.Entities.MeasurementCorrection> corrections)
		{
			if (corrections == null || corrections.Count == 0) return rawValue;

			if (rawValue < corrections[0].Measurement)
			{
				return rawValue + corrections[0].Correction;
			}

			int count = corrections.Count;
			for (int i = 1; i < count; i++)
			{
				if (rawValue < corrections[i].Measurement)
				{
					double d1 = rawValue - corrections[i-1].Measurement;
					double d2 = corrections[i].Measurement - rawValue;

					double corr;
					if (d1 + d2 <= float.Epsilon)
					{
						corr = (corrections[i - 1].Correction + corrections[i].Correction) / 2.0;
					}
					else
					{
						corr = (corrections[i - 1].Correction * d2 + corrections[i].Correction * d1) / (d1 + d2);
					}
					return rawValue + corr;
				}
			}

			return rawValue + corrections[count - 1].Correction;
		}


        /// <summary>
        /// Converts a measurement error to a correction (used when preparing correction tables).
        /// </summary>
        /// <param name="measuredValue">Measured value (in arbitrary units)</param>
        /// <param name="error">Measurement error in %</param>
        /// <returns>Correction in the same units as the measured value</returns>
        public static double CorrectionFromError(double measuredValue, double error)
        {
            double trueValue = measuredValue / (1 + error/100);
            double correction = trueValue - measuredValue;
            return correction;
        }


        /// <summary>
        /// Calculates the heat coefficient for water
        /// </summary>
        /// <param name="pressure">Pressure [bar]</param>
		/// <param name="T_in">Inlet temperature [°C]</param>
		/// <param name="T_out">Outlet temperature [°C]</param>
        /// <param name="flowMeasuredAtInlet">true = flow measured @inlet, false = flow measured @outlet</param>
        /// <returns> Heat coefficient for water [J/(m3 K)]</returns>
        public static double HeatCoefficientWater(double pressure, double T_in, double T_out, bool flowMeasuredAtInlet)
        {
            if (T_in == T_out) return 0;

            const double R = 461.526;                       /// [J kg^-1 K^-1]
            const double p_star_Pa = 16.53E6;               /// [Pa] (=16.53 MPa)
            const double T_star = 1386.0;                   /// [K]

            double T_in_K = Config.Units.ConvertTo(Config.Unit.K, T_in);
            double T_out_K = Config.Units.ConvertTo(Config.Unit.K, T_out);
            double tau_in = T_star / T_in_K;
            double tau_out = T_star / T_out_K;
			double pi = Config.Units.ConvertTo(Config.Unit.Pa, pressure) / p_star_Pa;

            double h_in = tau_in * GammaTau(pi, tau_in) * R * T_in_K;
            double h_out = tau_out * GammaTau(pi, tau_out) * R * T_out_K;

            double ni = flowMeasuredAtInlet ? GammaPi(pi, tau_in) * R * T_in_K / p_star_Pa
                                            : GammaPi(pi, tau_out) * R * T_out_K / p_star_Pa;

            return (h_in - h_out) / (ni * (T_in - T_out));
        }

        /// <summary>
        /// gamma(pi) see also STN EN 1434-1 Annex A (A.4)
        /// </summary>
        /// <param name="pi">pi = p / p* where p* = 16.53 MPa</param>
        /// <param name="tau">tau = T* / T where T* = 1386 K</param>
        /// <returns>gamma(pi)</returns>
        static double GammaPi(double pi, double tau)
        {
            double result = 0;
            for (int i = 0; i < 34; i++)
            {
                result -= ni[i] * Ii[i] * Math.Pow(7.1 - pi, Ii[i] - 1) * Math.Pow(tau - 1.222, Ji[i]);
            }
            return result;
        }

        /// <summary>
        /// gamma(tau) see also STN EN 1434-1 Annex A (A.7)
        /// </summary>
        /// <param name="pi">pi = p / p* where p* = 16.53 MPa</param>
        /// <param name="tau">tau = T* / T where T* = 1386 K</param>
        /// <returns>gamma(tau)</returns>
        static double GammaTau(double pi, double tau)
        {
            double result = 0;
            for (int i = 0; i < 34; i++)
            {
                result += ni[i] * Math.Pow(7.1 - pi, Ii[i]) * Ji[i] * Math.Pow(tau - 1.222, Ji[i] - 1);
            }
            return result;
        }

		/// <summary>
		/// Conversion of measured resistance of a platinum thermometer to temperature according to ITS-90
		/// </summary>
		/// <param name="R">Measured resistance in [°C]</param>
		/// <param name="R001C">Calibrated resistance in Ohm at 0.01°C</param>
		/// <param name="a7">Calibrated ITS-90 coefficient a7</param>
		/// <param name="b7">Calibrated ITS-90 coefficient b7</param>
		/// <param name="c7">Calibrated ITS-90 coefficient c7</param>
        /// <returns>Temperature in [°C]</returns>
		public static double PlatinumResistanceTM_ITS90_R2T(double R, double R001C, double a7, double b7, double c7)
		{
			double w = R / R001C;			/// ratio
			double r1 = w - 1.0;
			double dw = r1 * (a7 + r1 * (b7 + r1 * c7));	/// = a7*r1 + b7*r1^2 + c7*r1^3
			double wr = w - dw;
			double x = (wr - 2.64) / 1.64;

			double sum = 0;
			for (int i = Di.Length - 1; i >= 0; i--)
			{
				sum = sum * x + Di[i];
			}

			return sum;
		}

        /// <summary>
        /// Conversion of measured resistance of a platinum thermometer to temperature using Callendar-Van Dusen equations
        /// </summary>
        /// <param name="R">Measured resistance in [°C]</param>
        /// <param name="R0">Calibrated resistance in Ohm at 0°C</param>
        /// <param name="A">Calibration coefficient a</param>
        /// <param name="B">Calibration coefficient b</param>
        /// <returns>Temperature in [°C]</returns>
        public static double PlatinumResistanceTM_ITS27_R2T(double R, double R0, double A, double B)
        {
            if (R0 * R0 * A * A - 4 * R0 * B * (R0 - R) <= 0) return 0;     /// Out of range

            return (-(R0 * A) + Math.Sqrt(R0 * R0 * A * A - 4 * R0 * B * (R0 - R))) / (2 * R0 * B);
        }


		public static double Buoyancy()
		{
			return 1.00103;
		}
	}
}