Attempting to generate sound, and set cells automatically

Audio/RingBuffer.cs

@@ -0,0 +1,41 @@
+namespace FluidSim.Audio
+{
+    internal class RingBuffer
+    {
+        private readonly float[] buffer;
+        private volatile int readPos;
+        private volatile int writePos;
+
+        public RingBuffer(int capacity)
+        {
+            if ((capacity & (capacity - 1)) != 0)
+                throw new ArgumentException("Capacity must be a power of two.");
+            buffer = new float[capacity];
+        }
+
+        public int Count => (writePos - readPos) & (buffer.Length - 1);
+        public int Space => (readPos - writePos - 1) & (buffer.Length - 1);
+
+        public int Write(float[] data, int count)
+        {
+            int space = Space;
+            int toWrite = Math.Min(count, space);
+            int mask = buffer.Length - 1;
+            for (int i = 0; i < toWrite; i++)
+                buffer[(writePos + i) & mask] = data[i];
+            writePos = (writePos + toWrite) & mask;
+            return toWrite;
+        }
+
+        public int Read(float[] destination, int count)
+        {
+            int available = Count;
+            int toRead = Math.Min(count, available);
+            int mask = buffer.Length - 1;
+            for (int i = 0; i < toRead; i++)
+                destination[i] = buffer[(readPos + i) & mask];
+            readPos = (readPos + toRead) & mask;
+            return toRead;
+        }
+    }
+}

Audio/RingBufferStream.cs

@@ -0,0 +1,43 @@
+using SFML.Audio;
+using SFML.System;
+
+namespace FluidSim.Audio
+{
+    internal class RingBufferStream : SoundStream
+    {
+        private readonly RingBuffer ringBuffer;
+
+        public RingBufferStream(RingBuffer buffer)
+        {
+            ringBuffer = buffer;
+            // 2 channels, 44.1 kHz, standard stereo mapping
+            Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
+        }
+
+        protected override bool OnGetData(out short[] samples)
+        {
+            const int monoBlockSize = 512;          // number of mono samples we'll read
+            float[] temp = new float[monoBlockSize];
+            int read = ringBuffer.Read(temp, monoBlockSize);
+            samples = new short[monoBlockSize * 2];
+
+            if (read > 0)
+            {
+                for (int i = 0; i < read; i++)
+                {
+                    float clamped = Math.Clamp(temp[i], -1f, 1f);
+                    short final = (short)(clamped * short.MaxValue);
+                    samples[i * 2] = final;   // left
+                    samples[i * 2 + 1] = final;   // right
+                }
+            }
+            for (int i = read * 2; i < samples.Length; i++)
+                samples[i] = 0;
+
+            return true;
+        }
+
+        protected override void OnSeek(Time timeOffset) =>
+            throw new NotSupportedException();
+    }
+}

Audio/SoundEngine.cs

@@ -0,0 +1,45 @@
+namespace FluidSim.Audio;
+
+public class SoundEngine : IDisposable
+{
+    private readonly RingBuffer ringBuffer;
+    private readonly RingBufferStream stream;
+    private bool isPlaying;
+
+    public SoundEngine(int bufferCapacity = 16384)
+    {
+        ringBuffer = new RingBuffer(bufferCapacity);
+        stream = new RingBufferStream(ringBuffer);
+    }
+
+    public void Start()
+    {
+        if (isPlaying) return;
+        stream.Play();
+        isPlaying = true;
+    }
+
+    public void Stop()
+    {
+        if (!isPlaying) return;
+        stream.Stop();
+        isPlaying = false;
+        float[] drain = new float[ringBuffer.Count];
+        ringBuffer.Read(drain, drain.Length);
+    }
+
+    public int WriteSamples(float[] data, int count) =>
+        ringBuffer.Write(data, count);
+
+    public float Volume
+    {
+        get => stream.Volume;
+        set => stream.Volume = value;
+    }
+
+    public void Dispose()
+    {
+        Stop();
+        stream.Dispose();
+    }
+}

Audio/SoundProcessor.cs

@@ -0,0 +1,29 @@
+using FluidSim.Components;
+using FluidSim.Interfaces;
+using System;
+
+namespace FluidSim.Audio
+{
+    public static class SoundProcessor
+    {
+        public static float MaxDeltaP { get; set; } = 100_000f;
+        public static float MaxArea { get; set; } = 1e-4f;
+        public static float MaxVelocity { get; set; } = 343f;
+        public static float ReferenceDensity { get; set; } = 1.225f;
+        public static float ReferenceSpeedOfSound { get; set; } = 343f;
+        public static float Gain { get; set; } = 1.0f;
+
+        public static double ComputeSample(Connection conn)
+        {
+            Port portA = conn.PortA;
+            Port portB = conn.PortB;
+
+            double pressureUp = portA.Pressure;
+            double pressureDown = portB.Pressure;
+
+            // No flow or no pressure drop → silence
+            double deltaP = pressureUp - pressureDown;
+            return deltaP / 1;
+        }
+    }
+}

Components/Connection.cs

@@ -1,4 +1,6 @@
-namespace FluidSim.Interfaces
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
 {
     /// <summary>Pure data link between two ports, with orifice parameters.</summary>
     public class Connection

Components/Pipe1D.cs

@@ -3,85 +3,65 @@ using FluidSim.Interfaces;
 
 namespace FluidSim.Components
 {
-    public enum BoundaryType
-    {
-        VolumeCoupling,
-        OpenEnd,
-        ClosedEnd
-    }
-
     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;
-        private double _dx, _dt, _gamma, _area, _diameter;
+        private double _dx, _dt, _gamma = 1.4, _area;
         private double[] _rho, _rhou, _E;
+        private double _hydraulicDiameter;
 
-        // Volume‑coupling ghost states for boundaries A and B
+        private double _rhoLeft, _pLeft, _rhoRight, _pRight;
-        private double _rhoA, _pA;
+        private bool _leftBCSet, _rightBCSet;
-        private double _rhoB, _pB;
-        private bool _aBCSet, _bBCSet;
 
-        private BoundaryType _aBCType = BoundaryType.VolumeCoupling;
+        public double FrictionFactor { get; set; }
-        private BoundaryType _bBCType = BoundaryType.VolumeCoupling;
-
-        private double _aAmbientPressure = 101325.0;
-        private double _bAmbientPressure = 101325.0;
-
-        private const double CflTarget = 0.8;
-        private const double ReferenceSoundSpeed = 340.0;
 
         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);
 
-        public BoundaryType ABCType => _aBCType;
+        /// <summary>
-        public BoundaryType BBCType => _bBCType;
+        /// Create a pipe with CFL‑stable automatic cell count.
-
+        /// </summary>
-        public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
+        /// <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)
         {
-            double dtGlobal = 1.0 / sampleRate;
+            if (area <= 0) throw new ArgumentException("Pipe area must be > 0");
-            int nCells;
-
-            if (forcedCellCount > 1)
-            {
-                nCells = forcedCellCount;
-            }
-            else
-            {
-                double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
-                nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
-                while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
-                    nCells++;
-            }
-
-            _n = nCells;
-            _dx = length / _n;
-            _dt = dtGlobal;
             _area = area;
-            _gamma = 1.4;
+            _dt = 1.0 / sampleRate;
+            FrictionFactor = frictionFactor;
 
-            // Hydraulic diameter for a circular pipe
+            // Nyquist‑based cell count (wave resolution)
-            _diameter = 2.0 * Math.Sqrt(area / Math.PI);
+            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 SetABoundaryType(BoundaryType type) => _aBCType = type;
-        public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
-        public void SetAAmbientPressure(double p) => _aAmbientPressure = p;
-        public void SetBAmbientPressure(double p) => _bAmbientPressure = p;
-
         public void SetUniformState(double rho, double u, double p)
         {
             double e = p / ((_gamma - 1) * rho);
@@ -94,181 +74,23 @@ namespace FluidSim.Components
             }
         }
 
-        public void SetCellState(int i, double rho, double u, double p)
+        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)
         {
-            if (i < 0 || i >= _n) return;
+            _rhoLeft = rhoVol;
-            _rho[i] = rho;
+            _pLeft = pVol;
-            _rhou[i] = rho * u;
+            _leftBCSet = true;
-            double e = p / ((_gamma - 1) * rho);
-            _E[i] = rho * e + 0.5 * rho * u * u;
         }
 
-        public void SetAVolumeState(double rhoVol, double pVol)
+        public void SetRightVolumeState(double rhoVol, double pVol)
         {
-            _rhoA = rhoVol;
+            _rhoRight = rhoVol;
-            _pA = pVol;
+            _pRight = pVol;
-            _aBCSet = true;
+            _rightBCSet = true;
-        }
-
-        public void SetBVolumeState(double rhoVol, double pVol)
-        {
-            _rhoB = rhoVol;
-            _pB = pVol;
-            _bBCSet = true;
-        }
-
-        public void ClearBC() => _aBCSet = _bBCSet = false;
-
-        public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
-        {
-            double maxW = 0.0;
-            for (int i = 0; i < _n; i++)
-            {
-                double rho = _rho[i];
-                double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
-                double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
-                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)));
-        }
-
-        public void SimulateSingleStep(double dtSub)
-        {
-            int n = _n;
-            double[] Fm = new double[n + 1];
-            double[] Fp = new double[n + 1];
-            double[] Fe = new double[n + 1];
-
-            // ---------- Boundary A (face 0, left) ----------
-            double rhoIntA = _rho[0];
-            double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12);
-            double pIntA = Pressure(0);
-
-            switch (_aBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_aBCSet)
-                    {
-                        HLLCFlux(_rhoA, 0.0, _pA,
-                                 rhoIntA, uIntA, pIntA,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    }
-                    else
-                    {
-                        Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false,
-                                  out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-            }
-
-            // ---------- Internal faces ----------
-            for (int i = 0; i < n - 1; i++)
-            {
-                double rhoL = _rho[i];
-                double uL = _rhou[i] / Math.Max(rhoL, 1e-12);
-                double pL = Pressure(i);
-
-                double rhoR = _rho[i + 1];
-                double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12);
-                double pR = Pressure(i + 1);
-
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR,
-                         out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
-            }
-
-            // ---------- Boundary B (face n, right) ----------
-            double rhoIntB = _rho[n - 1];
-            double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12);
-            double pIntB = Pressure(n - 1);
-
-            switch (_bBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_bBCSet)
-                    {
-                        HLLCFlux(rhoIntB, uIntB, pIntB,
-                                 _rhoB, 0.0, _pB,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    }
-                    else
-                    {
-                        Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true,
-                                  out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-            }
-
-            // ---- Cell update with linear laminar damping ----
-            double radius = _diameter / 2.0;
-            double mu_air = 1.8e-5;
-            double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius);
-
-            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] -= dtSub * dM;
-                _rhou[i] -= dtSub * dP;
-                _E[i] -= dtSub * dE;
-
-                double rho = Math.Max(_rho[i], 1e-12);
-                double dampingRate = laminarCoeff / rho;
-                double dampingFactor = Math.Exp(-dampingRate * dtSub);
-                _rhou[i] *= dampingFactor;
-
-                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
-                double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
-                double pMin = 100.0;
-                double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
-                if (_E[i] < eMin) _E[i] = eMin;
-            }
-
-            // ---------- Port quantities ----------
-            double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0;
-            double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0;
-
-            PortA.MassFlowRate = mdotA_sub;
-            PortB.MassFlowRate = mdotB_sub;
-            PortA.Pressure = pIntA;
-            PortB.Pressure = pIntB;
-            PortA.Density = _rho[0];
-            PortB.Density = _rho[n - 1];
-
-            // Corrected enthalpy for both directions
-            if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet)
-            {
-                PortA.SpecificEnthalpy = mdotA_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(0)
-                    : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12);
-            }
-            if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet)
-            {
-                PortB.SpecificEnthalpy = mdotB_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(_n - 1)
-                    : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12);
-            }
         }
 
         private double GetCellTotalSpecificEnthalpy(int i)
@@ -280,114 +102,125 @@ namespace FluidSim.Components
             return h + 0.5 * u * u;
         }
 
-        private double Pressure(int i) =>
+        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));
 
-        // ========== Characteristic‑based Open End ==========
+        void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
-        private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb,
+                      out double fm, out double fp, out double fe)
-                                  out double fm, out double fp, out double fe)
         {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
+            const double eps = 1e-12;
+            pL = Math.Max(pL, eps);
+            pR = Math.Max(pR, eps);
 
-            // Subsonic inflow (uInt ≤ 0, so flow inside pipe ←)
+            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, eps));
-            if (uInt <= -cInt)   // supersonic inflow – use interior state as ghost
+            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, eps));
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt <= 0)   // subsonic inflow
-            {
-                // Reservoir condition: p = pAmb, T = 300K, u = 0
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-                return;
-            }
-            else   // subsonic outflow (uInt > 0)
-            {
-                // Ghost pressure forced to pAmb
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary)
-                double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0);
-
-                // Prevent spurious inflow by clipping to zero
-                if (uGhost < 0) uGhost = 0;
-
-                HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-            }
-        }
-
-        private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb,
-                                  out double fm, out double fp, out double fe)
-        {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
-
-            if (uInt >= cInt)   // supersonic outflow (extrapolation)
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt >= 0)   // subsonic outflow
-            {
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary)
-                double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0);
-
-                // Clip to zero to prevent inflow
-                if (uGhost > 0) uGhost = 0;
-
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe);
-            }
-            else   // subsonic inflow
-            {
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe);
-            }
-        }
-
-        // ========== Closed end (mirror) ==========
-        private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary,
-                                   out double fm, out double fp, out double fe)
-        {
-            double rhoGhost = rhoInt;
-            double pGhost = pInt;
-            double uGhost = -uInt;   // mirror velocity
-
-            if (isRightBoundary)
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
-            else
-                HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
-        }
-
-        // ========== Standard HLLC flux ==========
-        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 / Math.Max(rL, 1e-12));
-            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12));
             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 Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR))
+            double denom = rL * (SL - uL) - rR * (SR - uR);
-                      / (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;
@@ -396,25 +229,22 @@ namespace FluidSim.Components
             else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
             else if (Ss >= 0)
             {
-                double rsL = rL * (SL - uL) / (SL - Ss);
+                double diffSL = SL - uL;
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
+                if (Math.Abs(diffSL) < eps) diffSL = eps;
-                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
+                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 rsR = rR * (SR - uR) / (SR - Ss);
+                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 * (SR - uR)));
+                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * diffSR));
                 fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
             }
         }
-
-        public double GetPressureAtFraction(double fraction)
-        {
-            int i = (int)(fraction * (_n - 1));
-            i = Math.Clamp(i, 0, _n - 1);
-            return Pressure(i);
-        }
     }
 }

Components/SoundConnection.cs

@@ -0,0 +1,9 @@
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+    public class SoundConnection : Connection
+    {
+        public SoundConnection(Port a, Port b) : base(a, b) { }
+    }
+}

Components/Volume0D.cs

@@ -1,6 +1,4 @@
-using System;
+using FluidSim.Interfaces;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
 
 namespace FluidSim.Components
 {
@@ -50,26 +48,13 @@ namespace FluidSim.Components
             Port.SpecificEnthalpy = SpecificEnthalpy;
         }
 
-        // Original integrate (uses the constructor’s sample rate)
         public void Integrate()
-        {
-            Integrate(_dt);
-        }
-
-        public void SetPressure(double newPressure)
-        {
-            InternalEnergy = newPressure * Volume / (Gamma - 1.0);
-            // Mass stays the same, so density is unchanged
-        }
-
-        // New overload: integrate with a custom time step (for sub‑steps)
-        public void Integrate(double dtOverride)
         {
             double mdot = Port.MassFlowRate;
             double h_in = Port.SpecificEnthalpy;
 
-            double dm = mdot * dtOverride;
+            double dm = mdot * _dt;
-            double dE = (mdot * h_in) * dtOverride - Pressure * dVdt * dtOverride;
+            double dE = (mdot * h_in) * _dt - Pressure * dVdt * _dt;
 
             Mass += dm;
             InternalEnergy += dE;

Core/OrificeBoundary.cs

@@ -1,5 +1,4 @@
-using System;
+using FluidSim.Components;
-using FluidSim.Interfaces;
 
 namespace FluidSim.Core
 {

Core/Solver.cs

@@ -1,5 +1,4 @@
-using System;
+using FluidSim.Audio;
-using System.Collections.Generic;
 using FluidSim.Components;
 using FluidSim.Interfaces;
 
@@ -11,160 +10,116 @@ namespace FluidSim.Core
         private readonly List<Pipe1D> _pipes = new();
         private readonly List<Connection> _connections = new();
 
-        private double _dt;
+        public float LastSample { get; private set; }
 
         public void AddVolume(Volume0D v) => _volumes.Add(v);
         public void AddPipe(Pipe1D p) => _pipes.Add(p);
         public void AddConnection(Connection c) => _connections.Add(c);
-        public void SetTimeStep(double dt) => _dt = dt;
 
-        /// <summary>
+        public void Step()
-        /// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right).
-        /// </summary>
-        public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0)
         {
-            if (isA)
+            // 1. Publish volume states to their own ports
-            {
-                pipe.SetABoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetAAmbientPressure(ambientPressure);
-            }
-            else
-            {
-                pipe.SetBBoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetBAmbientPressure(ambientPressure);
-            }
-        }
-
-        public float Step()
-        {
-            // 1. Volumes publish state
             foreach (var v in _volumes)
                 v.PushStateToPort();
 
-            // 2. Set volume BCs for volume‑coupled ends
+            // 2. Handle direct volume‑to‑volume connections
+            foreach (var conn in _connections)
+            {
+                if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB))
+                {
+                    Volume0D volA = _volumes.Find(v => v.Port == conn.PortA);
+                    Volume0D volB = _volumes.Find(v => v.Port == conn.PortB);
+
+                    if (volA == null || volB == null) continue;
+
+                    double pA = volA.Pressure, rhoA = volA.Density;
+                    double pB = volB.Pressure, rhoB = volB.Density;
+
+                    double mdot = OrificeBoundary.MassFlow(pA, rhoA, pB, rhoB, conn);
+
+                    if (mdot > 0) // A → B
+                    {
+                        volA.Port.MassFlowRate = -mdot;
+                        volB.Port.MassFlowRate = mdot;
+                        volB.Port.SpecificEnthalpy = volA.SpecificEnthalpy; // fluid from A
+                        volA.Port.SpecificEnthalpy = volA.SpecificEnthalpy; // outflow carries its own enthalpy
+                    }
+                    else // B → A (mdot negative)
+                    {
+                        volA.Port.MassFlowRate = -mdot;  // positive
+                        volB.Port.MassFlowRate = mdot;  // negative
+                        volA.Port.SpecificEnthalpy = volB.SpecificEnthalpy; // fluid from B
+                        volB.Port.SpecificEnthalpy = volB.SpecificEnthalpy; // outflow carries its own
+                    }
+                }
+            }
+
+            // 3. Pipe‑volume boundary conditions – unchanged
             foreach (var conn in _connections)
             {
                 if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                {
+                    SetVolumeBC(conn.PortA, conn.PortB);
-                    var pipe = GetPipe(conn.PortA);
-                    bool isA = pipe.PortA == conn.PortA;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortA, conn.PortB);
-                }
                 else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                {
+                    SetVolumeBC(conn.PortB, conn.PortA);
-                    var pipe = GetPipe(conn.PortB);
-                    bool isA = pipe.PortB == conn.PortB;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortB, conn.PortA);
-                }
             }
 
-            // 3. Sub‑steps
+            // 4. Run pipe simulations
-            int nSub = 1;
             foreach (var p in _pipes)
-                nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
+                p.Simulate();
-            double dtSub = _dt / nSub;
 
-            for (int sub = 0; sub < nSub; sub++)
+            // 5. Transfer pipe‑to‑volume flows – unchanged
-            {
-                foreach (var p in _pipes)
-                    p.SimulateSingleStep(dtSub);
-
-                foreach (var conn in _connections)
-                {
-                    if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortA);
-                        bool isA = pipe.PortA == conn.PortA;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortA, conn.PortB, dtSub);
-                    }
-                    else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortB);
-                        bool isA = pipe.PortB == conn.PortB;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortB, conn.PortA, dtSub);
-                    }
-                }
-
-                if (sub < nSub - 1)
-                {
-                    foreach (var v in _volumes)
-                        v.PushStateToPort();
-
-                    foreach (var conn in _connections)
-                    {
-                        if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortA);
-                            bool isA = pipe.PortA == conn.PortA;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortA, conn.PortB);
-                        }
-                        else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortB);
-                            bool isA = pipe.PortB == conn.PortB;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortB, conn.PortA);
-                        }
-                    }
-                }
-            }
-
-            // 5. Audio samples (none for now, but placeholder)
-            var audioSamples = new List<float>();
             foreach (var conn in _connections)
             {
-                if (conn is SoundConnection sc)
+                if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                    audioSamples.Add(sc.GetAudioSample());
+                    TransferPipeToVolume(conn.PortA, conn.PortB);
+                else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+                    TransferPipeToVolume(conn.PortB, conn.PortA);
             }
 
-            // 6. Clear BC flags
+            // 6. Integrate volumes
-            foreach (var p in _pipes)
+            foreach (var v in _volumes)
-                p.ClearBC();
+                v.Integrate();
 
-            return SoundProcessor.MixAndClip(audioSamples.ToArray());
+            // 7. COMPUTE AUDIO SAMPLE from all sound connections
+            double sample = 0f;
+            foreach (var conn in _connections)
+            {
+                if (conn is SoundConnection)
+                {
+                    // Both ports have the same absolute mass flow; either works.
+                    sample += SoundProcessor.ComputeSample(conn);
+                }
+            }
+            LastSample = (float)Math.Tanh(sample);
         }
 
-        private bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
+        bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
-        private bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
+        bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
-        private Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
+        Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
-        private Volume0D GetVolume(Port p) => _volumes.Find(v => v.Port == p);
 
-        private void SetVolumeBC(Port pipePort, Port volPort)
+        void SetVolumeBC(Port pipePort, Port volPort)
         {
-            var pipe = GetPipe(pipePort);
+            Pipe1D pipe = GetPipe(pipePort);
             if (pipe == null) return;
-            bool isA = pipe.PortA == pipePort;
+            bool isLeft = pipe.PortA == pipePort;
-            if (isA)
+
-                pipe.SetAVolumeState(volPort.Density, volPort.Pressure);
+            if (isLeft)
+                pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure);
             else
-                pipe.SetBVolumeState(volPort.Density, volPort.Pressure);
+                pipe.SetRightVolumeState(volPort.Density, volPort.Pressure);
         }
 
-        private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub)
+        void TransferPipeToVolume(Port pipePort, Port volPort)
         {
             double mdot = pipePort.MassFlowRate;
             volPort.MassFlowRate = -mdot;
 
             if (mdot < 0)  // pipe → volume
             {
+                // pipePort.SpecificEnthalpy is already total (h + ½u²)
                 volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
             }
-            // else volume’s own enthalpy (from PushStateToPort) is used
+            // else: volume → pipe, volume’s own static enthalpy is used (already set)
-
-            GetVolume(volPort)?.Integrate(dtSub);
         }
     }
 }

Core/SoundEngine.cs

@@ -1,131 +0,0 @@
-using SFML.Audio;
-using SFML.System;
-
-namespace FluidSim;
-
-#region Lock‑free ring buffer (unchanged)
-internal class RingBuffer
-{
-    private readonly float[] buffer;
-    private volatile int readPos;
-    private volatile int writePos;
-
-    public RingBuffer(int capacity)
-    {
-        if ((capacity & (capacity - 1)) != 0)
-            throw new ArgumentException("Capacity must be a power of two.");
-        buffer = new float[capacity];
-    }
-
-    public int Count => (writePos - readPos) & (buffer.Length - 1);
-    public int Space => (readPos - writePos - 1) & (buffer.Length - 1);
-
-    public int Write(float[] data, int count)
-    {
-        int space = Space;
-        int toWrite = Math.Min(count, space);
-        int mask = buffer.Length - 1;
-        for (int i = 0; i < toWrite; i++)
-            buffer[(writePos + i) & mask] = data[i];
-        writePos = (writePos + toWrite) & mask;
-        return toWrite;
-    }
-
-    public int Read(float[] destination, int count)
-    {
-        int available = Count;
-        int toRead = Math.Min(count, available);
-        int mask = buffer.Length - 1;
-        for (int i = 0; i < toRead; i++)
-            destination[i] = buffer[(readPos + i) & mask];
-        readPos = (readPos + toRead) & mask;
-        return toRead;
-    }
-}
-#endregion
-
-#region Stereo stream that consumes the ring buffer
-internal class RingBufferStream : SoundStream
-{
-    private readonly RingBuffer ringBuffer;
-
-    public RingBufferStream(RingBuffer buffer)
-    {
-        ringBuffer = buffer;
-        // 2 channels, 44.1 kHz, standard stereo mapping
-        Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
-    }
-
-    protected override bool OnGetData(out short[] samples)
-    {
-        const int monoBlockSize = 512;          // number of mono samples we'll read
-        float[] temp = new float[monoBlockSize];
-        int read = ringBuffer.Read(temp, monoBlockSize);
-        samples = new short[monoBlockSize * 2];
-
-        if (read > 0)
-        {
-            for (int i = 0; i < read; i++)
-            {
-                float clamped = Math.Clamp(temp[i], -1f, 1f);
-                short final = (short)(clamped * short.MaxValue);
-                samples[i * 2] = final;   // left
-                samples[i * 2 + 1] = final;   // right
-            }
-        }
-        for (int i = read * 2; i < samples.Length; i++)
-            samples[i] = 0;
-
-        return true;
-    }
-
-    protected override void OnSeek(Time timeOffset) =>
-        throw new NotSupportedException();
-}
-#endregion
-
-#region Public sound engine API (unchanged)
-public class SoundEngine : IDisposable
-{
-    private readonly RingBuffer ringBuffer;
-    private readonly RingBufferStream stream;
-    private bool isPlaying;
-
-    public SoundEngine(int bufferCapacity = 16384)
-    {
-        ringBuffer = new RingBuffer(bufferCapacity);
-        stream = new RingBufferStream(ringBuffer);
-    }
-
-    public void Start()
-    {
-        if (isPlaying) return;
-        stream.Play();
-        isPlaying = true;
-    }
-
-    public void Stop()
-    {
-        if (!isPlaying) return;
-        stream.Stop();
-        isPlaying = false;
-        float[] drain = new float[ringBuffer.Count];
-        ringBuffer.Read(drain, drain.Length);
-    }
-
-    public int WriteSamples(float[] data, int count) =>
-        ringBuffer.Write(data, count);
-
-    public float Volume
-    {
-        get => stream.Volume;
-        set => stream.Volume = value;
-    }
-
-    public void Dispose()
-    {
-        Stop();
-        stream.Dispose();
-    }
-}
-#endregion

Core/SoundProcessor.cs

@@ -1,23 +0,0 @@
-namespace FluidSim.Core
-{
-    /// <summary>
-    /// Mixes multiple audio samples and applies a soft‑clipping tanh.
-    /// </summary>
-    public static class SoundProcessor
-    {
-        /// <summary>Overall gain applied after mixing (before tanh).</summary>
-        public static float MasterGain { get; set; } = 0.01f;
-
-        /// <summary>
-        /// Mixes an array of raw audio samples and returns a single sample in [‑1, 1].
-        /// </summary>
-        public static float MixAndClip(params float[] samples)
-        {
-            float sum = 0f;
-            foreach (float s in samples)
-                sum += s;
-            sum *= MasterGain;
-            return sum;
-        }
-    }
-}

Interfaces/SoundConnection.cs

@@ -1,25 +0,0 @@
-namespace FluidSim.Interfaces
-{
-    /// <summary>
-    /// A Connection that also produces an audio sample from the pressure drop across it.
-    /// </summary>
-    public class SoundConnection : Connection
-    {
-        /// <summary>Gain applied to the normalised pressure difference.</summary>
-        public float Gain { get; set; } = 1.0f;
-
-        /// <summary>Reference pressure used for normalisation (Pa). Default: 1 atm.</summary>
-        public double ReferencePressure { get; set; } = 101325.0;
-
-        public SoundConnection(Port a, Port b) : base(a, b) { }
-
-        /// <summary>
-        /// Returns a normalised audio sample proportional to the pressure difference.
-        /// </summary>
-        public float GetAudioSample()
-        {
-            double dp = PortA.Pressure - PortB.Pressure;
-            return (float)(dp / ReferencePressure) * Gain;
-        }
-    }
-}

Program.cs

@@ -2,183 +2,95 @@
 using SFML.Window;
 using SFML.System;
 using System.Diagnostics;
-using FluidSim.Core;
+using FluidSim.Scenarios;
+using FluidSim.Audio;
 
 namespace FluidSim;
 
 public class Program
 {
     private const int SampleRate = 44100;
-    private const double DrawFrequency = 60.0;
-    private static Scenario scenario;
-
-    // Speed control
-    //private static double desiredSpeed = 1.0;
-    private static double desiredSpeed = 0.0001;
-    private static double currentSpeed = desiredSpeed;
-    private const double MinSpeed = 0.0001;
-    private const double MaxSpeed = 1.0;
-    private const double ScrollFactor = 1.1;
-
-    // Space‑toggle state
-    private static double lastDesiredSpeed = 0.1;   // remembers the last non‑1.0 scroll speed
-    private static bool isRealTime = true;          // true when desiredSpeed == 1.0
-
     private static volatile bool running = true;
+    // Global step counter – incremented every simulation step
+    private static long stepCount = 0;
 
     public static void Main()
     {
         var mode = new VideoMode(new Vector2u(1280, 720));
-        var window = new RenderWindow(mode, "Pipe Resonator");
+        var window = new RenderWindow(mode, "Fluid Simulation");
         window.SetVerticalSyncEnabled(true);
         window.Closed += (_, _) => { running = false; window.Close(); };
-        window.MouseWheelScrolled += OnMouseWheel;
-        window.KeyPressed += OnKeyPressed;
 
-        var soundEngine = new SoundEngine(bufferCapacity: 16384);
+        var soundEngine = new SoundEngine(bufferCapacity: 2048);
         soundEngine.Volume = 70;
         soundEngine.Start();
 
-        //scenario = new PipeResonatorScenario();
-        //scenario = new HelmholtzResonatorScenario();
-        scenario = new SodShockTubeScenario();
-
-        scenario.Initialize(SampleRate);
-
         var stopwatch = Stopwatch.StartNew();
-        double lastDrawTime = 0.0;
-        double drawInterval = 1.0 / DrawFrequency;
-        double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds;
 
-        // Resampling buffer
+        // --- Warmup: fill audio buffer with silence ---
-        List<float> simBuffer = new List<float>(4096);
+        int warmupSamples = SampleRate / 2;                // 0.5 s
-        double readIndex = 0.0;
+        float[] warmup = new float[warmupSamples];
+        for (int i = 0; i < warmupSamples; i++)
+            warmup[i] = 0;
+        soundEngine.WriteSamples(warmup, warmupSamples);
 
-        for (int i = 0; i < 4; i++)
+        // Reset timer after warmup – this is the “real‑time zero”
-            simBuffer.Add(scenario.Process());
+        stopwatch.Restart();
+        stepCount = 0;                                     // simulation steps start now
 
-        long totalSimSteps = simBuffer.Count;
+        // --- Initialise the simulation scenario ---
-        long totalOutputSamples = 0;
+        Simulation.Initialize(SampleRate);
 
-        double lastRealTime = stopwatch.Elapsed.TotalSeconds;
+        const int chunkSize = 2048;
-        const int outputChunk = 256;
+        float[] buffer = new float[chunkSize];
-        float[] outputBuf = new float[outputChunk];
+
+        double lastLogTime = 0.0;                         // for periodic speed printout
 
         while (window.IsOpen)
         {
             window.DispatchEvents();
 
-            double currentRealTime = stopwatch.Elapsed.TotalSeconds;
+            // --- Compute how many audio samples are needed since last frame ---
-            double dtSpeed = currentRealTime - lastSpeedUpdateTime;
+            double currentTime = stopwatch.Elapsed.TotalSeconds;
-            lastSpeedUpdateTime = currentRealTime;
+            double elapsed = currentTime;                  // since stopwatch was reset
+            int samplesNeeded = (int)(elapsed * SampleRate) - (int)(stepCount);
+            // (stepCount is total generated samples, so we just need the remainder)
 
-            // Smoothly transition currentSpeed → desiredSpeed
+            // --- Generate the required number of simulation steps ---
-            // When toggling, desiredSpeed jumps, but currentSpeed follows with a smooth lerp
+            while (samplesNeeded > 0 && running)
-            double smoothingRate = 8.0;   // higher = faster catch‑up
-            currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed));
-
-            // ---------- Generate audio ----------
-            double targetAudioClock = currentRealTime + 0.05;
-
-            while (totalOutputSamples < targetAudioClock * SampleRate && running)
             {
-                int toGenerate = (int)Math.Min(
+                int toGenerate = Math.Min(samplesNeeded, chunkSize);
-                    (long)outputChunk,
-                    (long)(targetAudioClock * SampleRate) - totalOutputSamples
-                );
-                if (toGenerate <= 0) break;
-
-                double maxIndex = readIndex + (toGenerate - 1) * currentSpeed + 2;
-                int requiredSimIndex = (int)Math.Ceiling(maxIndex);
-                while (simBuffer.Count - 1 < requiredSimIndex)
-                {
-                    simBuffer.Add(scenario.Process());
-                    totalSimSteps++;
-                }
-
                 for (int i = 0; i < toGenerate; i++)
                 {
-                    int i0 = (int)readIndex;
+                    buffer[i] = Simulation.Process();
-                    int i1 = i0 + 1;
+                    stepCount++;
-                    double frac = readIndex - i0;
-
-                    float y0 = simBuffer[Math.Clamp(i0, 0, simBuffer.Count - 1)];
-                    float y1 = simBuffer[Math.Clamp(i1, 0, simBuffer.Count - 1)];
-                    outputBuf[i] = (float)(y0 + (y1 - y0) * frac);
-
-                    readIndex += currentSpeed;
-
-                    while (readIndex >= 1.0 && simBuffer.Count > 2)
-                    {
-                        simBuffer.RemoveAt(0);
-                        readIndex -= 1.0;
-                    }
                 }
-
+                soundEngine.WriteSamples(buffer, toGenerate);
-                int accepted = soundEngine.WriteSamples(outputBuf, toGenerate);
+                samplesNeeded -= toGenerate;
-                totalOutputSamples += accepted;
-
-                if (accepted < toGenerate)
-                    break;
             }
 
-            // ---------- Drawing & title ----------
+            // --- Display speed ---
-            if (currentRealTime - lastDrawTime >= drawInterval)
+            double simTime = stepCount / (double)SampleRate;
+            double wallTime = stopwatch.Elapsed.TotalSeconds;
+            double speed = (wallTime > 0) ? simTime / wallTime : 0.0;
+
+            // Update window title with instant speed
+            window.SetTitle($"FluidSim | Speed: {speed:F3}× | Steps: {stepCount}");
+
+            // Console log once per second
+            if (wallTime - lastLogTime >= 1.0)
             {
-                double actualSpeed = totalOutputSamples / (currentRealTime * SampleRate);
+                Console.WriteLine($"Speed: {speed:F3}×  ({stepCount} steps, {wallTime:F2}s wall)");
-                double simTime = totalSimSteps / (double)SampleRate;
+                lastLogTime = wallTime;
-                string toggleHint = isRealTime ? "[Space] slow mo" : "[Space] real time";
-                window.SetTitle(
-                    $"{toggleHint}  Sim: {simTime:F2}s | " +
-                    $"Speed: {currentSpeed:F4}x → {desiredSpeed:F4}x | " +
-                    $"Actual: {actualSpeed:F2}x"
-                );
-
-                window.Clear(Color.Black);
-                scenario.Draw(window);
-                window.Display();
-                lastDrawTime = currentRealTime;
             }
+
+            // --- Rendering (placeholder) ---
+            window.Clear(Color.Black);
+            window.Display();
         }
 
+        // --- Cleanup ---
         soundEngine.Dispose();
         window.Dispose();
     }
-
-    private static void OnMouseWheel(object? sender, MouseWheelScrollEventArgs e)
-    {
-        bool wasRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
-
-        if (e.Delta > 0)
-            desiredSpeed *= ScrollFactor;
-        else if (e.Delta < 0)
-            desiredSpeed /= ScrollFactor;
-
-        desiredSpeed = Math.Clamp(desiredSpeed, MinSpeed, MaxSpeed);
-
-        // Update the remembered slow-mo speed (unless we are exactly at 1.0)
-        if (!wasRealTime || Math.Abs(desiredSpeed - 1.0) > 1e-6)
-            lastDesiredSpeed = desiredSpeed;
-
-        // Update isRealTime flag
-        isRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
-    }
-
-    private static void OnKeyPressed(object? sender, KeyEventArgs e)
-    {
-        if (e.Code == Keyboard.Key.Space)
-        {
-            if (isRealTime)
-            {
-                // Switch to the remembered slow speed
-                desiredSpeed = lastDesiredSpeed;
-            }
-            else
-            {
-                // Switch back to real time
-                desiredSpeed = 1.0;
-            }
-            isRealTime = !isRealTime;
-        }
-    }
 }

Scenarios/HelmholtzResonatorScenario.cs

@@ -1,133 +0,0 @@
-using System;
-using FluidSim.Components;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
-using SFML.Graphics;
-using SFML.System;
-
-namespace FluidSim.Core
-{
-    public class HelmholtzResonatorScenario : Scenario
-    {
-        private Solver solver;
-        private Volume0D cavity;
-        private Pipe1D neck;
-        private Connection coupling;
-        private int stepCount;
-        private double time;
-        private double dt;
-        private double ambientPressure = 1.0 * Units.atm;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            // 1‑litre cavity, 10% over‑pressure
-            double cavityVolume = 1e-3;
-            double initialCavityPressure = 1.1 * ambientPressure;
-            cavity = new Volume0D(cavityVolume, initialCavityPressure, 300.0, sampleRate)
-            {
-                Gamma = 1.4,
-                GasConstant = 287.0
-            };
-
-            // Neck: length 10 cm, radius 1 cm
-            double neckLength = 0.1;
-            double neckRadius = 0.01;
-            double neckArea = Math.PI * neckRadius * neckRadius;
-            neck = new Pipe1D(neckLength, neckArea, sampleRate, forcedCellCount: 40);
-            neck.SetUniformState(1.225, 0.0, ambientPressure);
-
-            coupling = new Connection(neck.PortA, cavity.Port)
-            {
-                Area = neckArea,
-                DischargeCoefficient = 0.62,
-                Gamma = 1.4
-            };
-
-            solver = new Solver();
-            solver.SetTimeStep(dt);
-            solver.AddVolume(cavity);
-            solver.AddPipe(neck);
-            solver.AddConnection(coupling);
-
-            // Port A (left) = volume coupling, Port B (right) = open end
-            solver.SetPipeBoundary(neck, isA: true, BoundaryType.VolumeCoupling);
-            solver.SetPipeBoundary(neck, isA: false, BoundaryType.OpenEnd, ambientPressure);
-        }
-
-        public override float Process()
-        {
-            float sample = solver.Step();
-            time += dt;
-            stepCount++;
-
-            double pOpen = neck.GetCellPressure(neck.GetCellCount() - 1);
-            float audio = (float)((pOpen - ambientPressure) / ambientPressure);
-
-            if (stepCount % 20 == 0)
-            {
-                double pCav = cavity.Pressure;
-                double mdotA = neck.PortA.MassFlowRate;   // positive = into pipe (leaving cavity)
-                Console.WriteLine(
-                    $"t={time * 1e3:F2} ms step={stepCount} " +
-                    $"P_cav={pCav:F1} Pa, P_open={pOpen:F1} Pa, " +
-                    $"mdot_A={mdotA * 1e3:F4} g/s, audio={audio:F4}");
-            }
-
-            return audio;
-        }
-
-        public override void Draw(RenderWindow target)
-        {
-            float winW = target.GetView().Size.X;
-            float winH = target.GetView().Size.Y;
-            float centerY = winH / 2f;
-
-            // Cavity rectangle
-            float cavityWidth = 120f;
-            float cavityHeight = 180f;
-            var cavityRect = new RectangleShape(new Vector2f(cavityWidth, cavityHeight));
-            cavityRect.Position = new Vector2f(40f, centerY - cavityHeight / 2f);
-            cavityRect.FillColor = PressureColor(cavity.Pressure);
-            target.Draw(cavityRect);
-
-            // Neck drawn as tapered pipe
-            int n = neck.GetCellCount();
-            float neckStartX = 40f + cavityWidth + 10f;
-            float neckEndX = winW - 60f;
-            float neckLenPx = neckEndX - neckStartX;
-            float dx = neckLenPx / (n - 1);
-            float baseRadius = 20f;
-
-            Vertex[] vertices = new Vertex[n * 2];
-            for (int i = 0; i < n; i++)
-            {
-                float x = neckStartX + i * dx;
-                double p = neck.GetCellPressure(i);
-                float r = baseRadius * (float)(0.5 + 0.5 * Math.Tanh((p - ambientPressure) / (ambientPressure * 0.2)));
-                if (r < 4f) r = 4f;
-                Color col = PressureColor(p);
-                vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);
-                vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col);
-            }
-            target.Draw(vertices, PrimitiveType.TriangleStrip);
-
-            // Open end indicator
-            var arrow = new CircleShape(8f);
-            arrow.Position = new Vector2f(neckEndX - 4f, centerY - 4f);
-            arrow.FillColor = Color.White;
-            target.Draw(arrow);
-        }
-
-        private Color PressureColor(double pressure)
-        {
-            double range = ambientPressure * 0.1;
-            double t = Math.Clamp((pressure - ambientPressure) / range, -1.0, 1.0);
-            byte r = (byte)(t > 0 ? 255 * t : 0);
-            byte b = (byte)(t < 0 ? -255 * t : 0);
-            byte g = (byte)(255 * (1 - Math.Abs(t)));
-            return new Color(r, g, b);
-        }
-    }
-}

Scenarios/PipeResonatorScenario.cs

@@ -1,184 +0,0 @@
-using FluidSim.Components;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
-using SFML.Graphics;
-using SFML.System;
-using System;
-
-namespace FluidSim.Core
-{
-    public class PipeResonatorScenario : Scenario
-    {
-        private Solver solver;
-        private Pipe1D pipe;
-        private int stepCount;
-        private double time;
-        private double dt;
-        private double ambientPressure = 1.0 * Units.atm;
-        private bool enableLogging = true;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            double length = 2;
-            double radius = 50 * Units.mm;
-            double area = Units.AreaFromDiameter(radius);
-
-            pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: 80);
-            pipe.SetUniformState(1.225, 0.0, ambientPressure);
-
-            solver = new Solver();
-            solver.SetTimeStep(dt);
-            solver.AddPipe(pipe);
-            // Open end at port A (left), closed end at port B (right)
-            solver.SetPipeBoundary(pipe, isA: true, BoundaryType.OpenEnd, ambientPressure);
-            solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd);
-
-            // Initial pressure pulse
-            int pulseCells = 5;
-            double pulsePressure = 2 * ambientPressure;
-            for (int i = 0; i < pulseCells; i++)
-                pipe.SetCellState(i, 1.225, 0.0, pulsePressure);
-        }
-
-        public override float Process()
-        {
-            float sample = solver.Step();
-            time += dt;
-            stepCount++;
-
-            double pMid = pipe.GetPressureAtFraction(0.5);
-            sample = (float)((pMid - ambientPressure) / ambientPressure);
-
-            Log(sample);
-            return sample;
-        }
-
-        private void Log(float sample)
-        {
-            if (!enableLogging) return;
-            if (stepCount % 10 == 0 && stepCount < 1000)
-            {
-                double pMid = pipe.GetPressureAtFraction(0.5);
-                double pOpen = pipe.GetCellPressure(0);
-                double pClosed = pipe.GetCellPressure(pipe.GetCellCount() - 1);
-                Console.WriteLine(
-                    $"t = {time * 1e3:F3} ms  Step {stepCount:D4}: " +
-                    $"sample = {sample:F3}, " +
-                    $"P_mid = {pMid:F2} Pa ({pMid / ambientPressure:F4} atm), " +
-                    $"P_open = {pOpen:F2} Pa, P_closed = {pClosed:F2} Pa");
-            }
-        }
-
-        public override void Draw(RenderWindow target)
-        {
-            float winWidth = target.GetView().Size.X;
-            float winHeight = target.GetView().Size.Y;
-
-            float pipeCenterY = winHeight / 2f;
-            float margin = 60f;
-            float pipeStartX = margin;
-            float pipeEndX = winWidth - margin;
-            float pipeLengthPx = pipeEndX - pipeStartX;
-            int n = pipe.GetCellCount();
-            float dx = pipeLengthPx / (n - 1);   // spacing between cell centres
-
-            float baseRadius = 25f;
-            float rangeFactor = 1f;
-            float scaleFactor = 5f;
-
-            // ----- smoothstep helper -----
-            static float SmoothStep(float edge0, float edge1, float x)
-            {
-                float t = Math.Clamp((x - edge0) / (edge1 - edge0), 0f, 1f);
-                return t * t * (3f - 2f * t);
-            }
-
-            // ----- Pre‑compute cell positions and radii -----
-            var centers = new float[n];
-            var radii = new float[n];
-            for (int i = 0; i < n; i++)
-            {
-                double p = pipe.GetCellPressure(i);
-                float deviation = (float)Math.Tanh((p - ambientPressure) / ambientPressure / rangeFactor);
-                radii[i] = baseRadius * (1f + deviation * scaleFactor);
-                if (radii[i] < 2f) radii[i] = 2f;
-                centers[i] = pipeStartX + i * dx;
-            }
-
-            // ----- Build triangle‑strip vertices -----
-            int segmentsPerCell = 8;   // smoothness
-            int totalPoints = n + (n - 1) * segmentsPerCell;
-            Vertex[] stripVertices = new Vertex[totalPoints * 2];   // top + bottom for each point
-            int idx = 0;
-
-            for (int i = 0; i < n; i++)
-            {
-                // ---- Cell centre ----
-                float x = centers[i];
-                float r = radii[i];
-                double p = pipe.GetCellPressure(i);
-                Color col = PressureColor(p);
-
-                stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col);
-                stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col);
-
-                // ---- Intermediate segments after this cell (if not last) ----
-                if (i < n - 1)
-                {
-                    for (int s = 1; s <= segmentsPerCell; s++)
-                    {
-                        float t = s / (float)segmentsPerCell;
-                        float st = SmoothStep(0f, 1f, t);
-                        float xi = centers[i] + (centers[i + 1] - centers[i]) * t;
-                        float ri = radii[i] + (radii[i + 1] - radii[i]) * st;
-                        double pi = pipe.GetCellPressure(i) * (1 - t) + pipe.GetCellPressure(i + 1) * t;
-                        Color coli = PressureColor(pi);
-
-                        stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY - ri), coli);
-                        stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY + ri), coli);
-                    }
-                }
-            }
-
-            // Draw the pipe as a triangle strip
-            var pipeMesh = new VertexArray(PrimitiveType.TriangleStrip, (uint)stripVertices.Length);
-            for (int i = 0; i < stripVertices.Length; i++)
-                pipeMesh[(uint)i] = stripVertices[i];
-            target.Draw(pipeMesh);
-
-            // ----- Closed end indicator (right) -----
-            float wallThickness = 8f;
-            var wall = new RectangleShape(new Vector2f(wallThickness, winHeight * 0.6f));
-            wall.Position = new Vector2f(pipeEndX, pipeCenterY - winHeight * 0.6f / 2f);
-            wall.FillColor = new Color(180, 180, 180);
-            target.Draw(wall);
-        }
-
-        /// <summary>Blue (low) → Green (ambient) → Red (high).</summary>
-        private Color PressureColor(double pressure)
-        {
-            double range = ambientPressure * 0.05;  // ±5% gives full colour swing
-            double t = (pressure - ambientPressure) / range;
-            t = Math.Clamp(t, -1.0, 1.0);
-
-            byte r, g, b;
-            if (t < 0)
-            {
-                double factor = -t;
-                r = 0;
-                g = (byte)(255 * (1 - factor));
-                b = (byte)(255 * factor);
-            }
-            else
-            {
-                double factor = t;
-                r = (byte)(255 * factor);
-                g = (byte)(255 * (1 - factor));
-                b = 0;
-            }
-            return new Color(r, g, b);
-        }
-    }
-}

Scenarios/Scenario.cs

@@ -1,23 +0,0 @@
-using SFML.Graphics;
-
-namespace FluidSim.Core
-{
-    public abstract class Scenario
-    {
-        /// <summary>
-        /// Initialize the scenario with a given audio sample rate.
-        /// </summary>
-        public abstract void Initialize(int sampleRate);
-
-        /// <summary>
-        /// Advance one simulation step and return an audio sample.
-        /// The step size is 1 / sampleRate seconds.
-        /// </summary>
-        public abstract float Process();
-
-        /// <summary>
-        /// Draw the current simulation state onto the given SFML render target.
-        /// </summary>
-        public abstract void Draw(RenderWindow target);
-    }
-}

Scenarios/Simulation.cs

@@ -0,0 +1,93 @@
+using FluidSim.Audio;
+using FluidSim.Components;
+using FluidSim.Core;
+using FluidSim.Utilities;
+
+namespace FluidSim.Scenarios
+{
+    public static class Simulation
+    {
+        private static Solver solver;
+        private static Volume0D cavity, ambient;
+        private static Pipe1D neck;
+        private static Connection connNeckCavity, connNeckAmbient;
+        private static int stepCount;
+        private static double time;
+        private static double dt;
+
+        public static void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            // --- Cavity (the “bottle”) ---
+            double V_cav = 1.0 * Units.L;                 // 1 litre
+            cavity = new Volume0D(V_cav, 10 * Units.atm, Units.Celsius(20), sampleRate);
+
+            // --- Ambient (a huge “room”) ---
+            double V_amb = 1000.0;                        // 1000 m³ ≈ constant pressure
+            ambient = new Volume0D(V_amb, 1.0 * Units.atm, Units.Celsius(20), sampleRate);
+
+            // --- Neck (short pipe) ---
+            double length_neck = 80.0 * Units.mm;
+            double diam_neck = 10 * Units.mm;
+            double area_neck = Units.AreaFromRadius(diam_neck, Units.mm);                             // few cells enough for a short neck
+            neck = new Pipe1D(length_neck, area_neck, sampleRate);
+            neck.SetUniformState(ambient.Density, 0.0, ambient.Pressure);
+            neck.FrictionFactor = 0.02;                    // slight damping
+
+            // --- Connections ---
+            // Cavity-to-neck (full area)
+            connNeckCavity = new Connection(cavity.Port, neck.PortA)
+            {
+                Area = area_neck,
+                DischargeCoefficient = 1.0,
+                Gamma = 1.4
+            };
+            // Neck-to-ambient (SoundConnection to capture the radiated tone)
+            connNeckAmbient = new SoundConnection(neck.PortB, ambient.Port)
+            {
+                Area = area_neck,
+                DischargeCoefficient = 1.0,
+                Gamma = 1.4
+            };
+
+            // --- Solver ---
+            solver = new Solver();
+            solver.AddVolume(cavity);
+            solver.AddVolume(ambient);
+            solver.AddPipe(neck);
+            solver.AddConnection(connNeckCavity);
+            solver.AddConnection(connNeckAmbient);
+
+            // --- Sound tuning ---
+            SoundProcessor.MaxDeltaP = 0.1f * (float)Units.atm; // small Δp expected
+            SoundProcessor.MaxArea = (float)area_neck;
+            SoundProcessor.MaxVelocity = 343f;
+            SoundProcessor.ReferenceDensity = 1.2f;
+            SoundProcessor.ReferenceSpeedOfSound = 343f;
+            SoundProcessor.Gain = 10.0f;  // amplify because Δp is small
+        }
+
+        public static float Process()
+        {
+            solver.Step();
+            time += dt;
+            stepCount++;
+            Log();
+            return solver.LastSample;
+        }
+
+        public static void Log()
+        {
+            if (stepCount <= 200 || stepCount % (int)(0.5 / dt) == 0)
+            {
+                Console.WriteLine(
+                    $"t = {time * 1e3:F3} ms  " +
+                    $"Sample = {solver.LastSample:F6},  " +
+                    $"P_cav = {cavity.Pressure / 1e5:F6} bar,  " +
+                    $"flow_cav = {cavity.Port.MassFlowRate / 1e3:F4} g/s,  " +
+                    $"flow_neck = {neck.PortB.MassFlowRate * 1e3:F4} g/s");
+            }
+        }
+    }
+}

Scenarios/SodShockTubeScenario.cs

@@ -1,158 +0,0 @@
-using System;
-using FluidSim.Components;
-using FluidSim.Utils;
-using SFML.Graphics;
-using SFML.System;
-
-namespace FluidSim.Core
-{
-    public class SodShockTubeScenario : Scenario
-    {
-        private Solver solver;
-        private Pipe1D pipe;
-        private int stepCount;
-        private double time;
-        private double dt;
-        private double ambientPressure = 1.0 * Units.atm;
-        private const double GasConstant = 287.0;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-            double length = 1.0;
-            double area = 1.0;
-            int nCells = 200;
-
-            pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: nCells);
-            pipe.SetUniformState(0.125, 0.0, 0.1 * ambientPressure);  // right state
-
-            // Left half high pressure
-            for (int i = 0; i < nCells / 2; i++)
-                pipe.SetCellState(i, 1.0, 0.0, ambientPressure);
-
-            solver = new Solver();
-            solver.SetTimeStep(dt);
-            solver.AddPipe(pipe);
-            solver.SetPipeBoundary(pipe, isA: true, BoundaryType.ClosedEnd);
-            solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd);
-        }
-
-        public override float Process()
-        {
-            float sample = solver.Step();
-            time += dt;
-            stepCount++;
-
-            double pMid = pipe.GetPressureAtFraction(0.5);
-            float audio = (float)((pMid - ambientPressure) / ambientPressure);
-
-            bool log = true;
-
-            if (log)
-            {
-                int n = pipe.GetCellCount();
-                Console.WriteLine($"step {stepCount}:");
-                Console.WriteLine("i    rho (kg/m³)   p (Pa)       T (K)     u (m/s)");
-                for (int i = 0; i < n; i++)
-                {
-                    if (i % 10 == 0)
-                    {
-                        double rho = pipe.GetCellDensity(i);
-                        double p = pipe.GetCellPressure(i);
-                        double u = pipe.GetCellVelocity(i);
-                        double T = p / (rho * GasConstant);   // GasConstant = 287.0
-                        Console.WriteLine($"{i,-4} {rho,10:F4}   {p,10:F1}   {T,8:F2}   {u,10:F4}");
-                    }
-                }
-                Console.WriteLine();
-            }
-
-            return audio;
-        }
-
-        public override void Draw(RenderWindow target)
-        {
-            float winW = target.GetView().Size.X;
-            float winH = target.GetView().Size.Y;
-            float centerY = winH / 2f;
-            float margin = 40f;
-            float pipeStartX = margin;
-            float pipeEndX = winW - margin;
-            float pipeLenPx = pipeEndX - pipeStartX;
-            int n = pipe.GetCellCount();
-            float dx = pipeLenPx / (n - 1);
-            float baseRadius = 60f;
-
-            Vertex[] vertices = new Vertex[n * 2];
-            for (int i = 0; i < n; i++)
-            {
-                float x = pipeStartX + i * dx;
-
-                double p = pipe.GetCellPressure(i);
-                double rho = pipe.GetCellDensity(i);
-                double T = p / (rho * GasConstant);           // temperature in Kelvin
-
-                // Radius from pressure (exaggerated deviation)
-                float r = baseRadius * (float)(p / ambientPressure * 2);
-                if (r < 4f) r = 4f;
-
-                // Colour from temperature
-                Color col = TemperatureColor(T);
-
-                vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);
-                vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col);
-            }
-            target.Draw(vertices, PrimitiveType.TriangleStrip);
-
-            // Diaphragm marker (faint white line at initial interface)
-            float diaphragmX = pipeStartX + (n / 2) * dx;
-            var line = new RectangleShape(new Vector2f(2f, winH * 0.5f));
-            line.Position = new Vector2f(diaphragmX - 1f, centerY - winH * 0.25f);
-            line.FillColor = new Color(255, 255, 255, 80);
-            target.Draw(line);
-        }
-
-        /// <summary>
-        /// Custom temperature‑to‑hue mapping that matches the given Sod‑tube hue values:
-        ///   250 K → 176, 300 K → 122, 350 K → 120?, 450 K → 71.
-        /// Interpolates piecewise linearly, clamping outside [250,450].
-        /// </summary>
-        private Color TemperatureColor(double T)
-        {
-            // 1. Map temperature to hue (0–360)
-            double[] Tknots = { 250, 282, 353, 450 };
-            double[] Hknots = { 176, 179, 122, 71 };
-            double hue;
-            if (T <= Tknots[0]) hue = Hknots[0];
-            else if (T >= Tknots[^1]) hue = Hknots[^1];
-            else
-            {
-                int i = 0;
-                while (i < Tknots.Length - 1 && T > Tknots[i + 1]) i++;
-                double frac = (T - Tknots[i]) / (Tknots[i + 1] - Tknots[i]);
-                hue = Hknots[i] + frac * (Hknots[i + 1] - Hknots[i]);
-            }
-
-            // 2. Convert hue to RGB (S = 1, V = 1)
-            double h = hue / 60.0;
-            int sector = (int)Math.Floor(h);
-            double f = h - sector;
-            byte p = 0;
-            byte q = (byte)(255 * (1 - f));
-            byte tByte = (byte)(255 * f);
-            byte v = 255;
-
-            byte r, g, b;
-            switch (sector % 6)
-            {
-                case 0: r = v; g = tByte; b = p; break;
-                case 1: r = q; g = v; b = p; break;
-                case 2: r = p; g = v; b = tByte; break;
-                case 3: r = p; g = q; b = v; break;
-                case 4: r = tByte; g = p; b = v; break;
-                default: r = v; g = p; b = q; break;
-            }
-            return new Color(r, g, b);
-        }
-    }
-}

Sources/EffortSource.cs

@@ -1,10 +0,0 @@
-using System;
-using System.Collections.Generic;
-using System.Text;
-
-namespace FluidSim.Sources
-{
-    internal class EffortSource
-    {
-    }
-}

Sources/FlowSource.cs

@@ -1,10 +0,0 @@
-using System;
-using System.Collections.Generic;
-using System.Text;
-
-namespace FluidSim.Sources
-{
-    internal class FlowSource
-    {
-    }
-}

Utilities/Units.cs

@@ -1,6 +1,4 @@
-using System;
+namespace FluidSim.Utilities
-
-namespace FluidSim.Utils
 {
     public static class Units
     {
@@ -21,10 +19,10 @@ namespace FluidSim.Utils
 
         public static double Celsius(double tC) => tC + 273.15;
 
-        public static double AreaFromRadius(double radius) =>
+        public static double AreaFromRadius(double radius, double unit = mm) =>
-            Math.PI * (radius) * (radius);
+            Math.PI * (radius * unit) * (radius * unit);
 
-        public static double AreaFromDiameter(double diameter) =>
+        public static double AreaFromDiameter(double diameter, double unit = mm) =>
-            Math.PI * 0.25 * (diameter) * (diameter);
+            Math.PI * 0.25 * (diameter * unit) * (diameter * unit);
     }
 }