using System;
using FluidSim.Interfaces;

namespace FluidSim.Components
{
    public class Pipe1D
    {
        public Port PortA { get; }
        public Port PortB { get; }
        public double Area => _area;

        private int _n;
        private double _dx, _dt, _gamma = 1.4, _area;
        private double[] _rho, _rhou, _E;
        private double _hydraulicDiameter;

        private double _rhoLeft, _pLeft, _rhoRight, _pRight;
        private bool _leftBCSet, _rightBCSet;

        public double FrictionFactor { get; set; }

        public int GetCellCount() => _n;
        public double GetCellDensity(int i) => _rho[i];
        public double GetCellPressure(int i) => Pressure(i);
        public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);

        /// <summary>
        /// Create a pipe with CFL‑stable automatic cell count.
        /// </summary>
        /// <param name="length">Pipe length [m].</param>
        /// <param name="area">Cross‑sectional area [m²].</param>
        /// <param name="sampleRate">Simulation step rate [Hz].</param>
        /// <param name="c0">Speed of sound [m/s] (default 343).</param>
        /// <param name="frictionFactor">Darcy friction factor (0 = inviscid).</param>
        /// <param name="cflSafety">CFL safety factor ≤ 1 (0.8 recommended).</param>
        public Pipe1D(double length, double area, int sampleRate,
                      double c0 = 343.0, double frictionFactor = 0.02,
                      double cflSafety = 0.8)
        {
            if (area <= 0) throw new ArgumentException("Pipe area must be > 0");
            _area = area;
            _dt = 1.0 / sampleRate;
            FrictionFactor = frictionFactor;

            // Nyquist‑based cell count (wave resolution)
            double nNyquist = Math.Ceiling(length * sampleRate / c0);
            // CFL‑stable cell count: dx ≥ maxSpeed·dt / cflSafety, maxSpeed = 2·c0 (supersonic safe)
            double maxSpeed = 2.0 * c0;
            double dxMinStable = maxSpeed * _dt / cflSafety;
            double nStable = Math.Floor(length / dxMinStable);

            _n = Math.Max(2, (int)Math.Min(nNyquist, nStable));
            _dx = length / _n;

            _rho = new double[_n];
            _rhou = new double[_n];
            _E = new double[_n];

            _hydraulicDiameter = Math.Max(2.0 * Math.Sqrt(_area / Math.PI), 1e-9);

            PortA = new Port();
            PortB = new Port();
        }

        public void SetUniformState(double rho, double u, double p)
        {
            double e = p / ((_gamma - 1) * rho);
            double Etot = rho * e + 0.5 * rho * u * u;
            for (int i = 0; i < _n; i++)
            {
                _rho[i] = rho;
                _rhou[i] = rho * u;
                _E[i] = Etot;
            }
        }

        public double GetLeftPressure() => Pressure(0);
        public double GetRightPressure() => Pressure(_n - 1);
        public double GetLeftDensity() => _rho[0];
        public double GetRightDensity() => _rho[_n - 1];

        public void SetLeftVolumeState(double rhoVol, double pVol)
        {
            _rhoLeft = rhoVol;
            _pLeft = pVol;
            _leftBCSet = true;
        }

        public void SetRightVolumeState(double rhoVol, double pVol)
        {
            _rhoRight = rhoVol;
            _pRight = pVol;
            _rightBCSet = true;
        }

        private double GetCellTotalSpecificEnthalpy(int i)
        {
            double rho = Math.Max(_rho[i], 1e-12);
            double u = _rhou[i] / rho;
            double p = Pressure(i);
            double h = _gamma / (_gamma - 1.0) * p / rho;
            return h + 0.5 * u * u;
        }

        public void Simulate()
        {
            int n = _n;
            double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1];

            // --- Left boundary (face 0) ---
            if (_leftBCSet)
            {
                double rhoL = _rhoLeft, uL = 0.0, pL = _pLeft;
                double rhoR = _rho[0], uR = _rhou[0] / Math.Max(rhoR, 1e-12), pR = Pressure(0);
                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[0], out Fp[0], out Fe[0]);
            }
            else
            {
                Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0;
            }

            // --- Internal faces ---
            for (int i = 0; i < n - 1; i++)
            {
                double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
                double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
                HLLCFlux(_rho[i], uL, Pressure(i), _rho[i + 1], uR, Pressure(i + 1),
                         out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
            }

            // --- Right boundary (face n) ---
            if (_rightBCSet)
            {
                double rhoL = _rho[n - 1], uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12), pL = Pressure(n - 1);
                double rhoR = _rhoRight, uR = 0.0, pR = _pRight;
                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[n], out Fp[n], out Fe[n]);
            }
            else
            {
                Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0;
            }

            // --- Cell update (inviscid fluxes) ---
            for (int i = 0; i < n; i++)
            {
                double dM = (Fm[i + 1] - Fm[i]) / _dx;
                double dP = (Fp[i + 1] - Fp[i]) / _dx;
                double dE = (Fe[i + 1] - Fe[i]) / _dx;
                _rho[i] -= _dt * dM;
                _rhou[i] -= _dt * dP;
                _E[i] -= _dt * dE;

                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
                double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
                if (_E[i] < kinetic) _E[i] = kinetic;

                // Emergency reset if NaN
                if (double.IsNaN(_rho[i]) || double.IsNaN(_rhou[i]) || double.IsNaN(_E[i]))
                {
                    _rho[i] = 1.225;                   // reset to atmospheric air at 300 K
                    _rhou[i] = 0.0;
                    _E[i] = 101325.0 / (_gamma - 1.0); // internal energy at 1 atm
                }
            }

            // --- Friction (Darcy–Weisbach, energy‑conserving) ---
            if (FrictionFactor > 0)
            {
                double D = _hydraulicDiameter;
                double twoD = 2.0 * D;
                for (int i = 0; i < n; i++)
                {
                    double rho = _rho[i];
                    double u = _rhou[i] / rho;
                    double absU = Math.Abs(u);
                    double src = FrictionFactor * rho * absU * u / twoD;
                    double kinOld = 0.5 * rho * u * u;
                    _rhou[i] -= _dt * src;
                    double uNew = _rhou[i] / rho;
                    double kinNew = 0.5 * rho * uNew * uNew;
                    _E[i] += (kinOld - kinNew);
                }
            }

            // --- Publish to ports ---
            PortA.Pressure = Pressure(0);
            PortA.Density = _rho[0];
            PortB.Pressure = Pressure(_n - 1);
            PortB.Density = _rho[_n - 1];

            PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0;
            PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0;

            PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0);
            PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1);

            _leftBCSet = _rightBCSet = false;
        }

        double Pressure(int i) =>
            (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));

        void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
                      out double fm, out double fp, out double fe)
        {
            const double eps = 1e-12;
            pL = Math.Max(pL, eps);
            pR = Math.Max(pR, eps);

            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, eps));
            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, eps));
            double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL;
            double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR;

            double SL = Math.Min(uL - cL, uR - cR);
            double SR = Math.Max(uL + cL, uR + cR);

            double denom = rL * (SL - uL) - rR * (SR - uR);
            double Ss;
            if (Math.Abs(denom) < eps)
                Ss = 0.5 * (uL + uR);
            else
                Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;

            double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
            double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;

            if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
            else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
            else if (Ss >= 0)
            {
                double diffSL = SL - uL;
                if (Math.Abs(diffSL) < eps) diffSL = eps;
                double rsL = rL * diffSL / (SL - Ss);
                double ps = pL + rL * diffSL * (Ss - uL);
                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * diffSL));
                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
            }
            else
            {
                double diffSR = SR - uR;
                if (Math.Abs(diffSR) < eps) diffSR = eps;
                double rsR = rR * diffSR / (SR - Ss);
                double ps = pL + rL * (SL - uL) * (Ss - uL);
                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * diffSR));
                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
            }
        }
    }
}