using System;
using System.Numerics;
using FluidSim.Interfaces;

namespace FluidSim.Components
{
    public enum BoundaryType
    {
        OpenEnd,
        ZeroPressureOpen,
        ClosedEnd,
        GhostCell
    }

    public class Pipe1D
    {
        public Port PortA { get; }
        public Port PortB { get; }
        public double Area => _area;
        public double DampingMultiplier { get; set; } = 1.0;

        private int _n;                     // number of real cells
        private float _dx, _dt;            // spatial step, global time step
        private float _area, _diameter;    // cross‑section, diameter
        private float _gamma;              // ratio of specific heats (1.4)

        // Conserved variables – arrays sized [_n] (only real cells, ghosts handled externally)
        private float[] _rho;
        private float[] _rhou;
        private float[] _E;

        // Flux arrays for faces 0 .. _n (face i is between cell i-1 and i)
        private float[] _fluxM;   // mass flux
        private float[] _fluxP;   // momentum flux
        private float[] _fluxE;   // energy flux

        // Ghost cell states
        private float _rhoGhostL, _uGhostL, _pGhostL;
        private float _rhoGhostR, _uGhostR, _pGhostR;
        private bool _ghostLSet, _ghostRSet;

        private BoundaryType _aBCType = BoundaryType.GhostCell;
        private BoundaryType _bBCType = BoundaryType.GhostCell;

        private float _aAmbientPressure = 101325f;
        private float _bAmbientPressure = 101325f;

        // CFL / sub-stepping
        private const float CflTarget = 0.8f;
        private const float ReferenceSoundSpeed = 340f;
        private float _lastMassFlow = 0f;

        // Pre‑computed for damping
        private float _laminarCoeff;  // 8*mu / r^2, then multiplied by DampingMultiplier

        // ---- Energy loss (Newton cooling) ----
        private float _ambientEnergyReference;   // total energy density at ambient (Pamb / (γ-1))
        public float EnergyRelaxationRate { get; set; } = 0.0f;   // 1/s

        public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
        {
            float dtGlobal = 1f / sampleRate;
            int nCells;
            float dxTarget = ReferenceSoundSpeed * dtGlobal / CflTarget;

            if (forcedCellCount > 1)
                nCells = forcedCellCount;
            else
            {
                nCells = Math.Max(2, (int)Math.Round((float)length / dxTarget, MidpointRounding.AwayFromZero));
                while (length / nCells > dxTarget * 1.01f && nCells < int.MaxValue - 1)
                    nCells++;
            }

            _n = nCells;
            _dx = (float)(length / nCells);
            _dt = dtGlobal;
            _area = (float)area;
            _diameter = (float)(2.0 * Math.Sqrt(area / Math.PI));
            _gamma = 1.4f;

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

            // +1 because there are _n+1 faces
            _fluxM = new float[_n + 1];
            _fluxP = new float[_n + 1];
            _fluxE = new float[_n + 1];

            // Pre‑compute laminar damping coefficient (using air at 20°C)
            float mu_air = 1.8e-5f;
            float radius = _diameter * 0.5f;
            _laminarCoeff = 8f * mu_air / (radius * radius); // will be multiplied by DampingMultiplier at each step

            // Ambient reference energy (internal energy per unit volume at 101325 Pa)
            _ambientEnergyReference = 101325f / (_gamma - 1f);   // ≈ 253312.5 J/m³

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

        // ==================== PUBLIC API ============================
        public void SetABoundaryType(BoundaryType type) => _aBCType = type;
        public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
        public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p;
        public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p;

        public void SetGhostLeft(double rho, double u, double p)
        {
            _rhoGhostL = (float)rho;
            _uGhostL = (float)u;
            _pGhostL = (float)p;
            _ghostLSet = true;
        }
        public void SetGhostRight(double rho, double u, double p)
        {
            _rhoGhostR = (float)rho;
            _uGhostR = (float)u;
            _pGhostR = (float)p;
            _ghostRSet = true;
        }
        public void ClearGhostFlag()
        {
            _ghostLSet = false;
            _ghostRSet = false;
        }

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

        public int GetCellCount() => _n;
        public double GetCellDensity(int i) => _rho[i];
        public double GetCellVelocity(int i)
        {
            float rho = Math.Max(_rho[i], 1e-12f);
            return _rhou[i] / rho;
        }
        public double GetCellPressure(int i)
        {
            float rho = Math.Max(_rho[i], 1e-12f);
            return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
        }

        public double GetPressureAtFraction(double fraction)
        {
            int i = (int)(fraction * (_n - 1));
            i = Math.Clamp(i, 0, _n - 1);
            return GetCellPressure(i);
        }

        public void SetCellState(int i, double rho, double u, double p)
        {
            if (i < 0 || i >= _n) return;
            float r = (float)rho;
            float vel = (float)u;
            float pr = (float)p;
            _rho[i] = r;
            _rhou[i] = r * vel;
            float e = pr / ((_gamma - 1f) * r);
            _E[i] = r * e + 0.5f * r * vel * vel;
        }

        public double GetOpenEndMassFlow()
        {
            if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen)
                return 0.0;

            int lastCell = _n - 1;
            float rho = _rho[lastCell];
            float u = _rhou[lastCell] / Math.Max(rho, 1e-12f);
            float p = PressureScalar(lastCell);

            float c = MathF.Sqrt(_gamma * p / rho);
            float uFace = u;
            float rhoFace = rho;
            float pFace = p;

            if (uFace > 0 && uFace < c) // subsonic outflow
            {
                float s = p / MathF.Pow(rho, _gamma);
                float rhoAmb = MathF.Pow(_bAmbientPressure / s, 1f / _gamma);
                float cAmb = MathF.Sqrt(_gamma * _bAmbientPressure / rhoAmb);
                float J_plus = u + 2f * c / (_gamma - 1f);
                float uFaceNew = J_plus - 2f * cAmb / (_gamma - 1f);
                if (uFaceNew > 0) uFace = uFaceNew;
                if (uFace < 0) uFace = 0;
                rhoFace = rhoAmb;
                pFace = _bAmbientPressure;
            }

            return rhoFace * uFace * _area;
        }

        public double GetAndStoreMassFlowForDerivative()
        {
            float current = (float)GetOpenEndMassFlow();
            double derivative = (current - _lastMassFlow) / _dt;
            _lastMassFlow = current;
            return derivative;
        }

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

        // ==================== MAIN SIMULATION ==================================
        public void SimulateSingleStep(double dtSub)
        {
            float dt = (float)dtSub;
            int n = _n;

            // --- 1. Left boundary face (index 0) – scalar -----------------------
            float rhoL = _rho[0];
            float uL = _rhou[0] / Math.Max(rhoL, 1e-12f);
            float pL = PressureScalar(0);
            ComputeLeftBoundaryFlux(rhoL, uL, pL, out _fluxM[0], out _fluxP[0], out _fluxE[0]);

            // --- 2. Internal faces (1 .. n-1) – SIMD ---------------------------
            int vectorSize = Vector<float>.Count;
            int lastSimdFace = n - vectorSize;   // highest face index that starts a full vector block
            for (int f = 1; f <= lastSimdFace; f += vectorSize)
            {
                SimdInternalFluxBlock(f, vectorSize);
            }
            // Scalar remainder for faces f .. n-1
            for (int f = Math.Max(1, lastSimdFace + 1); f < n; f++)
            {
                float rhoLi = _rho[f - 1];
                float uLi = _rhou[f - 1] / Math.Max(rhoLi, 1e-12f);
                float pLi = PressureScalar(f - 1);
                float rhoRi = _rho[f];
                float uRi = _rhou[f] / Math.Max(rhoRi, 1e-12f);
                float pRi = PressureScalar(f);
                HLLCFluxScalar(rhoLi, uLi, pLi, rhoRi, uRi, pRi,
                               out _fluxM[f], out _fluxP[f], out _fluxE[f]);
            }

            // --- 3. Right boundary face (index n) – scalar --------------------
            float rhoR = _rho[n - 1];
            float uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12f);
            float pR = PressureScalar(n - 1);
            ComputeRightBoundaryFlux(rhoR, uR, pR, out _fluxM[n], out _fluxP[n], out _fluxE[n]);

            // --- 4. Cell update + damping + energy loss – SIMD -----------------
            SimdCellUpdate(dt);
        }

        // ==================== PRIVATE SCALAR HELPERS ===========================
        private float PressureScalar(int i)
        {
            float rho = Math.Max(_rho[i], 1e-12f);
            return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
        }

        private void ComputeLeftBoundaryFlux(float rhoInt, float uInt, float pInt,
            out float fm, out float fp, out float fe)
        {
            if (_aBCType == BoundaryType.GhostCell && _ghostLSet)
                HLLCFluxScalar(_rhoGhostL, _uGhostL, _pGhostL, rhoInt, uInt, pInt, out fm, out fp, out fe);
            else if (_aBCType == BoundaryType.OpenEnd)
                OpenEndFluxLeft(rhoInt, uInt, pInt, _aAmbientPressure, out fm, out fp, out fe);
            else if (_aBCType == BoundaryType.ZeroPressureOpen)
            {
                float rhoFace = rhoInt;
                float uFace = uInt;
                float pFace = _aAmbientPressure;
                HLLCFluxScalar(rhoFace, uFace, pFace, rhoInt, uInt, pInt, out fm, out fp, out fe);
            }
            else if (_aBCType == BoundaryType.ClosedEnd)
                ClosedEndFlux(rhoInt, uInt, pInt, false, out fm, out fp, out fe);
            else
            { fm = 0; fp = pInt; fe = 0; }
        }

        private void ComputeRightBoundaryFlux(float rhoInt, float uInt, float pInt,
            out float fm, out float fp, out float fe)
        {
            if (_bBCType == BoundaryType.GhostCell && _ghostRSet)
                HLLCFluxScalar(rhoInt, uInt, pInt, _rhoGhostR, _uGhostR, _pGhostR, out fm, out fp, out fe);
            else if (_bBCType == BoundaryType.OpenEnd)
                OpenEndFluxRight(rhoInt, uInt, pInt, _bAmbientPressure, out fm, out fp, out fe);
            else if (_bBCType == BoundaryType.ZeroPressureOpen)
            {
                float rhoFace = rhoInt;
                float uFace = uInt;
                float pFace = _bAmbientPressure;
                HLLCFluxScalar(rhoInt, uInt, pInt, rhoFace, uFace, pFace, out fm, out fp, out fe);
            }
            else if (_bBCType == BoundaryType.ClosedEnd)
                ClosedEndFlux(rhoInt, uInt, pInt, true, out fm, out fp, out fe);
            else
            { fm = 0; fp = pInt; fe = 0; }
        }

        // ==================== SCALAR HLLC & BOUNDARY FLUX ======================
        private void HLLCFluxScalar(float rL, float uL, float pL, float rR, float uR, float pR,
                                    out float fm, out float fp, out float fe)
        {
            float cL = MathF.Sqrt(_gamma * pL / Math.Max(rL, 1e-12f));
            float cR = MathF.Sqrt(_gamma * pR / Math.Max(rR, 1e-12f));
            float EL = pL / ((_gamma - 1f) * rL) + 0.5f * uL * uL;
            float ER = pR / ((_gamma - 1f) * rR) + 0.5f * uR * uR;
            float SL = Math.Min(uL - cL, uR - cR);
            float SR = Math.Max(uL + cL, uR + cR);

            float denom = rL * (SL - uL) - rR * (SR - uR);
            float Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;

            float FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
            float 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)
            {
                float rsL = rL * (SL - uL) / (SL - Ss);
                float ps = pL + rL * (SL - uL) * (Ss - uL);
                float EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
            }
            else
            {
                float rsR = rR * (SR - uR) / (SR - Ss);
                float ps = pR + rR * (SR - uR) * (Ss - uR);
                float EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
            }
        }

        private void OpenEndFluxLeft(float rhoInt, float uInt, float pInt, float pAmb,
                                     out float fm, out float fp, out float fe)
        {
            float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
            if (uInt <= -cInt) // supersonic inflow
            {
                fm = rhoInt * uInt;
                fp = rhoInt * uInt * uInt + pInt;
                fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
                return;
            }
            if (uInt <= 0) // subsonic inflow
            {
                float T0 = 300f, R = 287f;
                float ghostRho = pAmb / (R * T0);
                HLLCFluxScalar(ghostRho, 0f, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
                return;
            }
            // subsonic outflow
            float s = pInt / MathF.Pow(rhoInt, _gamma);
            float ghostRho2 = MathF.Pow(pAmb / s, 1f / _gamma);
            float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho2);
            float J_minus = uInt - 2f * cInt / (_gamma - 1f);
            float uGhost = J_minus + 2f * cGhost / (_gamma - 1f);
            if (uGhost < 0) uGhost = 0;
            HLLCFluxScalar(ghostRho2, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
        }

        private void OpenEndFluxRight(float rhoInt, float uInt, float pInt, float pAmb,
                                      out float fm, out float fp, out float fe)
        {
            float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
            if (uInt >= cInt) // supersonic outflow
            {
                fm = rhoInt * uInt;
                fp = rhoInt * uInt * uInt + pInt;
                fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
                return;
            }
            if (uInt >= 0) // subsonic outflow
            {
                float s = pInt / MathF.Pow(rhoInt, _gamma);
                float ghostRho = MathF.Pow(pAmb / s, 1f / _gamma);
                float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho);
                float J_plus = uInt + 2f * cInt / (_gamma - 1f);
                float uGhost = J_plus - 2f * cGhost / (_gamma - 1f);
                if (uGhost > 0) uGhost = 0;
                HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho, uGhost, pAmb, out fm, out fp, out fe);
                return;
            }
            // subsonic inflow
            float T0 = 300f, R = 287f;
            float ghostRho2 = pAmb / (R * T0);
            HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho2, 0f, pAmb, out fm, out fp, out fe);
        }

        private void ClosedEndFlux(float rhoInt, float uInt, float pInt, bool isRight,
                                   out float fm, out float fp, out float fe)
        {
            float rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt;
            if (isRight)
                HLLCFluxScalar(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
            else
                HLLCFluxScalar(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
        }

        // ==================== SIMD INTERNAL FACE ROUTINE ========================
        private void SimdInternalFluxBlock(int startFace, int count)
        {
            int leftIdx = startFace - 1;
            int rightIdx = startFace;

            Vector<float> rL = new Vector<float>(_rho, leftIdx);
            Vector<float> ruL = new Vector<float>(_rhou, leftIdx);
            Vector<float> EL = new Vector<float>(_E, leftIdx);

            Vector<float> rR = new Vector<float>(_rho, rightIdx);
            Vector<float> ruR = new Vector<float>(_rhou, rightIdx);
            Vector<float> ER = new Vector<float>(_E, rightIdx);

            Vector<float> uL = ruL / rL;
            Vector<float> uR = ruR / rR;

            Vector<float> half = new Vector<float>(0.5f);
            Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
            Vector<float> gammaVec = new Vector<float>(_gamma);

            Vector<float> pL = gammaMinus1 * (EL - half * ruL * ruL / rL);
            Vector<float> pR = gammaMinus1 * (ER - half * ruR * ruR / rR);

            Vector<float> cL = Vector.SquareRoot(gammaVec * pL / rL);
            Vector<float> cR = Vector.SquareRoot(gammaVec * pR / rR);

            Vector<float> SL = Vector.Min(uL - cL, uR - cR);
            Vector<float> SR = Vector.Max(uL + cL, uR + cR);

            Vector<float> num = (pR - pL) + rL * uL * (SL - uL) - rR * uR * (SR - uR);
            Vector<float> den = rL * (SL - uL) - rR * (SR - uR);
            Vector<float> Ss = num / den;

            Vector<float> eL = EL / rL;
            Vector<float> eR = ER / rR;

            // Left flux
            Vector<float> Fm_L = ruL;
            Vector<float> Fp_L = ruL * uL + pL;
            Vector<float> Fe_L = (EL + pL) * uL;

            // Right flux
            Vector<float> Fm_R = ruR;
            Vector<float> Fp_R = ruR * uR + pR;
            Vector<float> Fe_R = (ER + pR) * uR;

            // Star‑left fluxes
            Vector<float> diffL = SL - uL;
            Vector<float> dL_den = SL - Ss;
            Vector<float> rsL = rL * diffL / dL_den;
            Vector<float> psSL = pL + rL * diffL * (Ss - uL);
            Vector<float> EsL = eL + (Ss - uL) * (Ss + pL / (rL * diffL));
            Vector<float> Fm_starL = rsL * Ss;
            Vector<float> Fp_starL = rsL * Ss * Ss + psSL;
            Vector<float> Fe_starL = (rsL * EsL + psSL) * Ss;

            // Star‑right fluxes
            Vector<float> diffR = SR - uR;
            Vector<float> dR_den = SR - Ss;
            Vector<float> rsR = rR * diffR / dR_den;
            Vector<float> psSR = pR + rR * diffR * (Ss - uR);
            Vector<float> EsR = eR + (Ss - uR) * (Ss + pR / (rR * diffR));
            Vector<float> Fm_starR = rsR * Ss;
            Vector<float> Fp_starR = rsR * Ss * Ss + psSR;
            Vector<float> Fe_starR = (rsR * EsR + psSR) * Ss;

            var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector<float>.Zero);
            var maskSRle0 = Vector.LessThanOrEqual(SR, Vector<float>.Zero);
            var maskMiddle = ~(maskSLge0 | maskSRle0);
            var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector<float>.Zero);
            var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector<float>.Zero);

            Vector<float> fm = Vector.ConditionalSelect(maskSLge0, Fm_L,
                              Vector.ConditionalSelect(maskSRle0, Fm_R,
                              Vector.ConditionalSelect(maskStarL, Fm_starL,
                              Vector.ConditionalSelect(maskStarR, Fm_starR, Vector<float>.Zero))));

            Vector<float> fp = Vector.ConditionalSelect(maskSLge0, Fp_L,
                              Vector.ConditionalSelect(maskSRle0, Fp_R,
                              Vector.ConditionalSelect(maskStarL, Fp_starL,
                              Vector.ConditionalSelect(maskStarR, Fp_starR, Vector<float>.Zero))));

            Vector<float> fe = Vector.ConditionalSelect(maskSLge0, Fe_L,
                              Vector.ConditionalSelect(maskSRle0, Fe_R,
                              Vector.ConditionalSelect(maskStarL, Fe_starL,
                              Vector.ConditionalSelect(maskStarR, Fe_starR, Vector<float>.Zero))));

            fm.CopyTo(_fluxM, startFace);
            fp.CopyTo(_fluxP, startFace);
            fe.CopyTo(_fluxE, startFace);
        }

        // ==================== SIMD CELL UPDATE + DAMPING + ENERGY LOSS =========
        private void SimdCellUpdate(float dt)
        {
            float dt_dx = dt / _dx;
            Vector<float> vDtDx = new Vector<float>(dt_dx);
            float coeff = _laminarCoeff * (float)DampingMultiplier;
            Vector<float> vCoeff = new Vector<float>(coeff);
            Vector<float> vDt = new Vector<float>(dt);

            int vectorSize = Vector<float>.Count;
            int n = _n;
            int lastSimdCell = n - vectorSize;

            // Pre‑defined constants used in clamping
            Vector<float> half = new Vector<float>(0.5f);
            Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
            Vector<float> ambientEnergyVec = new Vector<float>(_ambientEnergyReference);
            Vector<float> energyRelaxRateVec = new Vector<float>(EnergyRelaxationRate);

            for (int i = 0; i <= lastSimdCell; i += vectorSize)
            {
                // Load conserved
                Vector<float> r = new Vector<float>(_rho, i);
                Vector<float> ru = new Vector<float>(_rhou, i);
                Vector<float> E = new Vector<float>(_E, i);

                // Load fluxes at faces i (left) and i+1 (right)
                Vector<float> flxM_L = new Vector<float>(_fluxM, i);
                Vector<float> flxM_R = new Vector<float>(_fluxM, i + 1);
                Vector<float> flxP_L = new Vector<float>(_fluxP, i);
                Vector<float> flxP_R = new Vector<float>(_fluxP, i + 1);
                Vector<float> flxE_L = new Vector<float>(_fluxE, i);
                Vector<float> flxE_R = new Vector<float>(_fluxE, i + 1);

                // Update conserved: Q_new = Q - dt/dx * (flux_right - flux_left)
                Vector<float> newR = r - vDtDx * (flxM_R - flxM_L);
                Vector<float> newRu = ru - vDtDx * (flxP_R - flxP_L);
                Vector<float> newE = E - vDtDx * (flxE_R - flxE_L);

                // Damping
                Vector<float> dampingFactor = Vector.Exp(-vCoeff / r * vDt);
                newRu *= dampingFactor;

                // Energy loss (Newton cooling toward ambient)
                Vector<float> relaxFactor = Vector.Exp(-energyRelaxRateVec * vDt);
                newE = ambientEnergyVec + (newE - ambientEnergyVec) * relaxFactor;

                // Clamp density
                newR = Vector.Max(newR, new Vector<float>(1e-12f));

                // Enforce minimal pressure (100 Pa) -> minimal energy
                Vector<float> kinE = half * newRu * newRu / newR;
                Vector<float> eMin = new Vector<float>(100f) / gammaMinus1 + kinE;
                newE = Vector.Max(newE, eMin);

                newR.CopyTo(_rho, i);
                newRu.CopyTo(_rhou, i);
                newE.CopyTo(_E, i);
            }

            // Scalar remainder
            float relaxRate = EnergyRelaxationRate;
            for (int i = Math.Max(0, lastSimdCell + 1); i < n; i++)
            {
                float r = _rho[i];
                float ru = _rhou[i];
                float E = _E[i];

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

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

                // Damping
                float dampingFactor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt);
                newRu *= dampingFactor;

                // Energy loss
                float relaxFactor = MathF.Exp(-relaxRate * dt);
                newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;

                // Clamps
                if (newR < 1e-12f) newR = 1e-12f;
                float kin = 0.5f * newRu * newRu / newR;
                float emin = 100f / (_gamma - 1f) + kin;
                if (newE < emin) newE = emin;

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