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); } /// /// 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]. /// 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); } } }