using System;
using FluidSim.Interfaces;

namespace FluidSim.Components
{
    /// <summary>
    /// 1‑D compressible Euler pipe (finite‑volume, HLLC flux).
    /// Boundary conditions are set externally via SetGhostLeft/Right.
    /// Enforces that ghosts are always valid before stepping.
    /// Uses exponential damping and Newtonian energy relaxation.
    /// </summary>
    public class Pipe1D : IComponent
    {
        public Port PortA { get; }
        public Port PortB { get; }
        public double Area { get; }
        public double DampingMultiplier { get; set; } = 1.0;
        public double EnergyRelaxationRate { get; set; } = 0.0;   // 1/s

        // Ambient pressure for the energy relaxation term (default 101325 Pa)
        private double _ambientPressure = 101325.0;
        public double AmbientPressure
        {
            get => _ambientPressure;
            set
            {
                _ambientPressure = value;
                _ambientEnergyReference = value / (_gamma - 1.0);
            }
        }

        // Geometry
        private readonly int _n;                // number of real cells
        private readonly double _dx;            // cell size (m)
        private readonly double _diameter;      // m
        private readonly double _gamma = 1.4;

        // Conserved variables [0 .. _n-1]
        private double[] _rho;
        private double[] _rhou;
        private double[] _E;

        // Face fluxes [0 .. _n]
        private double[] _fluxM;
        private double[] _fluxP;
        private double[] _fluxE;

        // Ghost cells (set externally)
        private double _rhoGhostL, _uGhostL, _pGhostL;
        private double _rhoGhostR, _uGhostR, _pGhostR;
        private bool _ghostLValid, _ghostRValid;

        // Pre‑computed damping coefficient
        private double _laminarCoeff;
        private double _ambientEnergyReference;   // internal energy density at ambient pressure

        /// <summary>
        /// Initialise a pipe with a given cell count.
        /// </summary>
        /// <param name="length">Pipe length (m).</param>
        /// <param name="area">Cross‑sectional area (m²).</param>
        /// <param name="cellCount">Number of finite‑volume cells (≥ 4).</param>
        public Pipe1D(double length, double area, int cellCount)
        {
            if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4");

            _n = cellCount;
            _dx = length / _n;
            Area = area;
            _diameter = 2.0 * Math.Sqrt(area / Math.PI);

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

            _fluxM = new double[_n + 1];
            _fluxP = new double[_n + 1];
            _fluxE = new double[_n + 1];

            // Laminar damping coefficient for air at 20°C (multiplied by DampingMultiplier each step)
            double mu_air = 1.8e-5;
            double radius = _diameter * 0.5;
            _laminarCoeff = 8.0 * mu_air / (radius * radius);

            // Ambient energy reference (internal energy per unit volume at 101325 Pa)
            _ambientEnergyReference = 101325.0 / (_gamma - 1.0);

            PortA = new Port { Owner = this };
            PortB = new Port { Owner = this };

            // Initial state = still air at ambient conditions
            SetUniformState(1.225, 0.0, 101325.0);
        }

        IReadOnlyList<Port> IComponent.Ports => new[] { PortA, PortB };

        // No integration needed for pipes – state is advanced via sub‑steps
        public void UpdateState(double dt) { }

        // ---------- Ghost cell interface ----------
        public void SetGhostLeft(double rho, double u, double p)
        {
            _rhoGhostL = rho;
            _uGhostL = u;
            _pGhostL = p;
            _ghostLValid = true;
        }

        public void SetGhostRight(double rho, double u, double p)
        {
            _rhoGhostR = rho;
            _uGhostR = u;
            _pGhostR = p;
            _ghostRValid = true;
        }

        public void ClearGhostFlags()
        {
            _ghostLValid = false;
            _ghostRValid = false;
        }

        public (double rho, double u, double p) GetInteriorStateLeft()
        {
            double rho = Math.Max(_rho[0], 1e-12);
            double u = _rhou[0] / rho;
            double p = PressureScalar(0);
            return (rho, u, p);
        }

        public (double rho, double u, double p) GetInteriorStateRight()
        {
            double rho = Math.Max(_rho[_n - 1], 1e-12);
            double u = _rhou[_n - 1] / rho;
            double p = PressureScalar(_n - 1);
            return (rho, u, p);
        }

        public int CellCount => _n;

        public double GetCellDensity(int i) => _rho[i];
        public double GetCellVelocity(int i)
        {
            double rho = Math.Max(_rho[i], 1e-12);
            return _rhou[i] / rho;
        }
        public double GetCellPressure(int i) => PressureScalar(i);

        // ---------- Sub‑stepping ----------
        public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
        {
            double maxW = 0.0;
            for (int i = 0; i < _n; i++)
            {
                double rho = Math.Max(_rho[i], 1e-12);
                double u = Math.Abs(_rhou[i] / rho);
                double p = PressureScalar(i);
                double c = Math.Sqrt(_gamma * p / rho);
                double local = u + c;
                if (local > maxW) maxW = local;
            }
            maxW = Math.Max(maxW, 1e-8);
            return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
        }

        // ---------- Main simulation step (per sub‑step) ----------
        public void SimulateSingleStep(double dtSub)
        {
            // Enforce that both ends have been provided with ghost states
            if (!_ghostLValid || !_ghostRValid)
                throw new InvalidOperationException("Pipe boundary ghosts not set before SimulateSingleStep.");

            double dt = dtSub;
            int n = _n;

            // Left boundary face (index 0)
            HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, _rho[0], _rhou[0] / _rho[0], PressureScalar(0),
                     out _fluxM[0], out _fluxP[0], out _fluxE[0]);

            // Internal faces (1 .. n-1)
            for (int f = 1; f < n; f++)
            {
                double rhoL = Math.Max(_rho[f - 1], 1e-12);
                double uL = _rhou[f - 1] / rhoL;
                double pL = PressureScalar(f - 1);
                double rhoR = Math.Max(_rho[f], 1e-12);
                double uR = _rhou[f] / rhoR;
                double pR = PressureScalar(f);
                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out _fluxM[f], out _fluxP[f], out _fluxE[f]);
            }

            // Right boundary face (index n)
            HLLCFlux(_rho[_n - 1], _rhou[_n - 1] / _rho[_n - 1], PressureScalar(_n - 1),
                     _rhoGhostR, _uGhostR, _pGhostR,
                     out _fluxM[n], out _fluxP[n], out _fluxE[n]);

            // Cell update
            double dt_dx = dt / _dx;
            double coeff = _laminarCoeff * DampingMultiplier;
            double relaxRate = EnergyRelaxationRate;

            for (int i = 0; i < n; i++)
            {
                double r = _rho[i];
                double ru = _rhou[i];
                double E = _E[i];

                double dM = _fluxM[i + 1] - _fluxM[i];
                double dP = _fluxP[i + 1] - _fluxP[i];
                double dE_flux = _fluxE[i + 1] - _fluxE[i];

                double newR = r - dt_dx * dM;
                double newRu = ru - dt_dx * dP;
                double newE = E - dt_dx * dE_flux;

                // Wall friction damping (laminar)
                double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
                newRu *= dampingFactor;

                // Newtonian cooling toward ambient energy
                double relaxFactor = Math.Exp(-relaxRate * dt);
                newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;

                // Clamps – minimum density 1e-12, minimum pressure 100 Pa
                newR = Math.Max(newR, 1e-12);
                double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12);
                double eMin = 100.0 / (_gamma - 1.0) + kin;
                newE = Math.Max(newE, eMin);

                _rho[i] = newR;
                _rhou[i] = newRu;
                _E[i] = newE;
            }

            // Update port states to reflect the current interior state (for audio / monitoring)
            (double rhoA, double uA, double pA) = GetInteriorStateLeft();
            PortA.Pressure = pA;
            PortA.Density = rhoA;
            PortA.Temperature = pA / (rhoA * 287.0);
            PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pA / rhoA;

            (double rhoB, double uB, double pB) = GetInteriorStateRight();
            PortB.Pressure = pB;
            PortB.Density = rhoB;
            PortB.Temperature = pB / (rhoB * 287.0);
            PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pB / rhoB;
        }

        // ---------- Private helpers ----------
        private double PressureScalar(int i)
        {
            double rho = Math.Max(_rho[i], 1e-12);
            return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho);
        }

        /// <summary>
        /// HLLC approximate Riemann solver (Toro, 1997).
        /// Computes the numerical flux at a face given left and right states.
        /// </summary>
        private void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
                              out double fm, out double fp, out double fe)
        {
            double cL = Math.Sqrt(_gamma * pL / rL);
            double cR = Math.Sqrt(_gamma * pR / rR);
            double EL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;   // specific total energy
            double ER = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;

            // Wave speed estimates (Davis, 1988)
            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 = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;

            double Fm_L = rL * uL;
            double Fp_L = rL * uL * uL + pL;
            double Fe_L = (rL * EL + pL) * uL;

            double Fm_R = rR * uR;
            double Fp_R = rR * uR * uR + pR;
            double Fe_R = (rR * ER + pR) * uR;

            if (SL >= 0) { fm = Fm_L; fp = Fp_L; fe = Fe_L; }
            else if (SR <= 0) { fm = Fm_R; fp = Fp_R; fe = Fe_R; }
            else if (Ss >= 0)
            {
                double rsL = rL * (SL - uL) / (SL - Ss);
                double ps = pL + rL * (SL - uL) * (Ss - uL);
                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
                fm = rsL * Ss;
                fp = rsL * Ss * Ss + ps;
                fe = (rsL * EsL + ps) * Ss;
            }
            else
            {
                double rsR = rR * (SR - uR) / (SR - Ss);
                double ps = pR + rR * (SR - uR) * (Ss - uR);
                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
                fm = rsR * Ss;
                fp = rsR * Ss * Ss + ps;
                fe = (rsR * EsR + ps) * Ss;
            }
        }

        /// <summary>Initialise all cells to a uniform state (rho, u, p).</summary>
        public void SetUniformState(double rho, double u, double p)
        {
            double e = p / ((_gamma - 1.0) * rho);
            double E = rho * e + 0.5 * rho * u * u;
            for (int i = 0; i < _n; i++)
            {
                _rho[i] = rho;
                _rhou[i] = rho * u;
                _E[i] = E;
            }
        }
    }
}