FluidSim/Scenarios/PipeResonatorScenario.cs

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 = 0.5;
            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);
            solver.SetPipeBoundary(pipe, isLeft: true, BoundaryType.OpenEnd, ambientPressure);
            solver.SetPipeBoundary(pipe, isLeft: 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);
        }
    }
}