diff --git a/Core/OutdoorExhaustReverb.cs b/Core/OutdoorExhaustReverb.cs new file mode 100644 index 0000000..3712771 --- /dev/null +++ b/Core/OutdoorExhaustReverb.cs @@ -0,0 +1,170 @@ +using System; + +namespace FluidSim.Core +{ + public class OutdoorExhaustReverb + { + // ---- Geometry ---- + private const float GroundReflDelay = 0.008f; // 8 ms (≈1.3 m) + private const float WallRefl1Delay = 0.045f; // ≈15 m + private const float WallRefl2Delay = 0.080f; // ≈27 m + + private DelayLine groundRefl; + private DelayLine wallRefl1; + private DelayLine wallRefl2; + + // ---- FDN for late diffuse tail ---- + private const int FDN_CHANNELS = 8; // dense, realistic + private DelayLine[] fdnDelays; + private float[] fdnState; + private OrthonormalMixer mixer; // energy‑preserving mixing + private LowPassFilter[] channelFilters; // per‑channel air absorption + + public float DryMix { get; set; } = 1.0f; + public float EarlyMix { get; set; } = 0.5f; + public float TailMix { get; set; } = 0.9f; + public float Feedback { get; set; } = 0.75f; // safe range 0.7‑0.9 + public float DampingFreq { get; set; } = 6000f; // Hz, above which air absorbs strongly + public float MatrixCoeff { get; set; } = 0.5f; // (kept for compatibility, not used) + + public OutdoorExhaustReverb(int sampleRate) + { + // Early reflection lines + groundRefl = new DelayLine((int)(sampleRate * GroundReflDelay)); + wallRefl1 = new DelayLine((int)(sampleRate * WallRefl1Delay)); + wallRefl2 = new DelayLine((int)(sampleRate * WallRefl2Delay)); + + // FDN delays: prime numbers for dense modal density (70‑150 ms) + int[] baseLengths = { 3203, 4027, 5521, 7027, 8521, 10007, 11503, 13009 }; + fdnDelays = new DelayLine[FDN_CHANNELS]; + for (int i = 0; i < FDN_CHANNELS; i++) + { + int len = Math.Min(baseLengths[i], (int)(sampleRate * 0.25)); // max 250 ms + fdnDelays[i] = new DelayLine(len); + } + + fdnState = new float[FDN_CHANNELS]; + mixer = new OrthonormalMixer(FDN_CHANNELS); + + // Air absorption: a gentle first‑order low‑pass per channel + channelFilters = new LowPassFilter[FDN_CHANNELS]; + float initialCutoff = DampingFreq; + for (int i = 0; i < FDN_CHANNELS; i++) + channelFilters[i] = new LowPassFilter(sampleRate, initialCutoff); + } + + public float Process(float drySample) + { + // ---- Early reflections ---- + float g = groundRefl.ReadWrite(drySample * 0.8f); + float w1 = wallRefl1.ReadWrite(drySample * 0.5f); + float w2 = wallRefl2.ReadWrite(drySample * 0.4f); + float early = (g + w1 + w2) * EarlyMix; + + // ---- FDN diffuse tail ---- + // Read the delayed outputs (which were stored last iteration) + float[] delOut = new float[FDN_CHANNELS]; + for (int i = 0; i < FDN_CHANNELS; i++) + delOut[i] = fdnDelays[i].Read(); + + // Mix the delayed outputs with the orthonormal matrix -> scattered signals + mixer.Process(delOut, fdnState); // result written into fdnState + + // Add fresh input to all channels + for (int i = 0; i < FDN_CHANNELS; i++) + fdnState[i] = drySample * 0.15f + Feedback * fdnState[i]; + + // Air absorption: per‑channel one‑pole low‑pass + for (int i = 0; i < FDN_CHANNELS; i++) + fdnState[i] = channelFilters[i].Process(fdnState[i]); + + // Write the new states into the delay lines + for (int i = 0; i < FDN_CHANNELS; i++) + fdnDelays[i].Write(fdnState[i]); + + // The tail output is the sum of the delayed outputs *before* the loop + float tailSum = 0f; + for (int i = 0; i < FDN_CHANNELS; i++) + tailSum += delOut[i]; + float tail = tailSum * TailMix; + + // Final mix + return drySample * DryMix + early + tail; + } + + // ---------- Helper classes (same as before but with separate Read/Write) ---------- + private class DelayLine + { + private float[] buffer; + private int writePos; + public DelayLine(int length) + { + buffer = new float[Math.Max(length, 1)]; + writePos = 0; + } + // Separated Read/Write to avoid ringing with immediate feedback + public float Read() + { + return buffer[writePos]; + } + public void Write(float value) + { + buffer[writePos] = value; + writePos = (writePos + 1) % buffer.Length; + } + // Old combined method (not used in FDN, only for early reflections) + public float ReadWrite(float value) + { + float outVal = buffer[writePos]; + buffer[writePos] = value; + writePos = (writePos + 1) % buffer.Length; + return outVal; + } + } + + private class LowPassFilter + { + private float b0, a1; + private float y1; + private float sampleRate; + public LowPassFilter(int sampleRate, float cutoff) + { + this.sampleRate = sampleRate; + SetCutoff(cutoff); + } + public void SetCutoff(float cutoff) + { + float w = 2 * (float)Math.PI * cutoff / sampleRate; + float a0 = 1 + w; + b0 = w / a0; + a1 = (1 - w) / a0; // first‑order low‑pass + } + public float Process(float x) + { + float y = b0 * x - a1 * y1; + y1 = y; + return y; + } + } + + /// + /// Computes a fast orthonormal mixing matrix (like Hadamard, but energy‑preserving). + /// + private class OrthonormalMixer + { + private int size; + public OrthonormalMixer(int size) { this.size = size; } + + public void Process(float[] input, float[] output) + { + // Simple energy‑conserving “allpass” mixing: + // Use a Householder reflection: y = (2/n) * sum(x) * ones - x + float sum = 0; + for (int i = 0; i < size; i++) sum += input[i]; + float factor = 2.0f / size; + for (int i = 0; i < size; i++) + output[i] = factor * sum - input[i]; + } + } + } +} \ No newline at end of file diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs index 3fb99bd..9229e9f 100644 --- a/Core/SoundProcessor.cs +++ b/Core/SoundProcessor.cs @@ -1,132 +1,48 @@ using System; +using FluidSim.Interfaces; namespace FluidSim.Core { public class SoundProcessor { - private double dt; - private double pipeArea; - private double ambientPressure = 101325.0; + private readonly double dt; + private readonly double scaleFactor; // 1 / (4π r) and a user gain + private double prevMassFlowOut; - // Monopole source state - private double lastMassFlow = 0.0; + // Simple low‑pass for derivative smoothing (≈ 2‑3 ms) + private double smoothDMdt; + private readonly double alpha; - // Gains - private float masterGain = 0.0005f; - private float pressureGain = 0.12f; - private float turbulenceGain = 0.05f; - private float turbulence = 0.05f; + public float Gain { get; set; } = 1.0f; - private PinkNoiseGenerator pinkNoise; - - // Reverb (outdoor) - private float[] delayLine; - private int writeIndex; - private float feedback = 0.50f; - private float lowpassCoeff = 0.70f; - private float lastFeedbackSample = 0f; - - public SoundProcessor(int sampleRate, double pipeDiameterMeters, float reverbTimeMs = 200.0f) + public SoundProcessor(int sampleRate, double listenerDistanceMeters = 1.0) { dt = 1.0 / sampleRate; - pipeArea = Math.PI * Math.Pow(pipeDiameterMeters / 2.0, 2.0); + scaleFactor = 1.0 / (4.0 * Math.PI * listenerDistanceMeters); - int delaySamples = (int)(sampleRate * reverbTimeMs / 1000.0); - delayLine = new float[delaySamples]; - writeIndex = 0; - - pinkNoise = new PinkNoiseGenerator(); + // Smoothing time constant ~ 2 ms, blocks single‑sample spikes + double tau = 0.002; + alpha = Math.Exp(-dt / tau); } - public float MasterGain + public float Process(Port port) { - get => masterGain; - set => masterGain = value; - } - public float PressureGain - { - get => pressureGain; - set => pressureGain = value; - } - public float TurbulenceGain - { - get => turbulenceGain; - set => turbulenceGain = value; - } - public float Turbulence - { - get => turbulence; - set => turbulence = value; - } + // Outflow mass flow (positive = leaving pipe) + double flowOut = -port.MassFlowRate; - public void SetAmbientPressure(double p) => ambientPressure = p; - public void SetPipeDiameter(double diameterMeters) => - pipeArea = Math.PI * Math.Pow(diameterMeters / 2.0, 2.0); + // Derivative + double rawDerivative = (flowOut - prevMassFlowOut) / dt; + prevMassFlowOut = flowOut; - public float Process(float massFlow, float pipeEndPressure) - { - // 1. Monopole: d(mdot)/dt - double derivative = (massFlow - lastMassFlow) / dt; - lastMassFlow = massFlow; - float monopole = (float)(derivative * masterGain); + // Smooth the derivative to kill isolated spikes + smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative; - // 2. Pressure component - float pressureDiff = (float)((pipeEndPressure - ambientPressure) / ambientPressure) * pressureGain; + // Far‑field monopole pressure + double pressure = smoothDMdt * scaleFactor * Gain; - float mixed = monopole + pressureDiff; - - // 3. Turbulence: amplitude ∝ U^8 - double velocity = massFlow / (pipeArea * 1.225); - double turbulenceAmp = Math.Pow(Math.Abs(velocity) * turbulence, 3.0); - float pink = pinkNoise.Next() * turbulenceGain * (float)turbulenceAmp; - - float combined = mixed + pink; - - // 4. Outdoor reverb - float delayed = delayLine[writeIndex]; - float filteredDelay = delayed * lowpassCoeff + lastFeedbackSample * (1f - lowpassCoeff); - lastFeedbackSample = filteredDelay; - float wet = delayed + filteredDelay * feedback; - delayLine[writeIndex] = combined + filteredDelay * feedback; - writeIndex = (writeIndex + 1) % delayLine.Length; - - // 5. Dry/wet mix - float output = combined * 0.7f + wet * 0.3f; - return MathF.Tanh(output); - } - } - - // PinkNoiseGenerator unchanged, same as before - internal class PinkNoiseGenerator - { - private readonly Random random = new Random(); - private readonly float[] whiteNoise = new float[7]; - private int currentIndex = 0; - - public PinkNoiseGenerator() - { - for (int i = 0; i < 7; i++) - whiteNoise[i] = (float)(random.NextDouble() * 2.0 - 1.0); - } - - public float Next() - { - whiteNoise[0] = (float)(random.NextDouble() * 2.0 - 1.0); - currentIndex = (currentIndex + 1) & 0x7F; - int updateMask = 0; - int temp = currentIndex; - for (int i = 0; i < 7; i++) - { - if ((temp & 1) == 0) - updateMask |= (1 << i); - temp >>= 1; - } - for (int i = 1; i < 7; i++) - if ((updateMask & (1 << i)) != 0) - whiteNoise[i] = (float)(random.NextDouble() * 2.0 - 1.0); - float sum = 0f; - for (int i = 0; i < 7; i++) sum += whiteNoise[i]; - return sum / 3.5f; + // Soft clip to ±1 for audio output (safe limit) + float sample = (float)Math.Tanh(pressure); + return sample; } } } \ No newline at end of file diff --git a/Scenarios/EngineScenario.cs b/Scenarios/EngineScenario.cs index 769ac2e..ceda9d2 100644 --- a/Scenarios/EngineScenario.cs +++ b/Scenarios/EngineScenario.cs @@ -1,5 +1,7 @@ using System; using FluidSim.Components; +using FluidSim.Utils; +using FluidSim.Interfaces; using SFML.Graphics; using SFML.System; @@ -13,108 +15,137 @@ namespace FluidSim.Core private Pipe1D exhaustPipe; private PipeVolumeConnection coupling; private SoundProcessor soundProcessor; + private OutdoorExhaustReverb reverb; + + private Port exitPort = new Port(); private double dt; - private double ambientPressure = 101325.0; + private double pipeArea; + private const double AmbientPressure = 101325.0; private double time; private int stepCount = 0; private const int LogInterval = 10000; - // Throttle 0..1 → target combustion pressure - public double Throttle { get; set; } = 0.05; // tiny throttle to keep idle - private const double IdlePeakPressure = 5.0 * 101325.0; // 5 bar - private const double MaxPeakPressure = 50.0 * 101325.0; // 50 bar + // Throttle 0..1 + public double Throttle { get; set; } = 0.0; // start with a light idle throttle + + // ---- Realistic combustion parameters ---- + private const double FullLoadPeakPressure = 70.0 * 101325.0; // 15 bar + + // ---- Idle speed governor ---- + private const double TargetIdleRPM = 800.0; // rad/s = RPM * π/30, we'll convert public override void Initialize(int sampleRate) { dt = 1.0 / sampleRate; - // Crankshaft (inertia + friction) - crankshaft = new Crankshaft(initialRPM: 100.0) // starter speed + // ---- Crankshaft: inertia + friction that gives ~800 RPM at idle ---- + crankshaft = new Crankshaft(initialRPM: 600.0) // start a bit low { - Inertia = 0.05, - FrictionConstant = 1.0, - FrictionViscous = 0.01 + Inertia = 0.005, // slightly heavier flywheel + FrictionConstant = 0.8, // static friction + FrictionViscous = 0.01 // viscous (linear with RPM) }; - // Pipe - double pipeLength = 0.5; + // ---- Pipe: add a tiny bit of damping to prevent unrealistic shocks ---- + double pipeLength = 2; double pipeRadius = 0.1; - double pipeArea = Math.PI * pipeRadius * pipeRadius; + pipeArea = Math.PI * pipeRadius * pipeRadius; exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 60); - exhaustPipe.SetUniformState(1.225, 0.0, ambientPressure); - exhaustPipe.EnergyRelaxationRate = 0f; - exhaustPipe.DampingMultiplier = 0; + exhaustPipe.SetUniformState(1.225, 0.0, AmbientPressure); + exhaustPipe.DampingMultiplier = 5; + exhaustPipe.EnergyRelaxationRate = 50; - // Cylinder (coupled to crankshaft) + // ---- Cylinder ---- engineCyl = new EngineCylinder(crankshaft, bore: 0.065, stroke: 0.0565, compressionRatio: 10.0, pipeArea: pipeArea, sampleRate: sampleRate); - // Coupling (valve → pipe) + // ---- Coupling ---- coupling = new PipeVolumeConnection(engineCyl.Cylinder, exhaustPipe, true, orificeArea: 0.0); - // Solver + // ---- Solver ---- solver = new Solver(); solver.SetTimeStep(dt); solver.AddVolume(engineCyl.Cylinder); solver.AddPipe(exhaustPipe); solver.AddConnection(coupling); - solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, ambientPressure); + solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, AmbientPressure); - // Sound (your tuned gains) - soundProcessor = new SoundProcessor(sampleRate, pipeRadius * 2, reverbTimeMs: 500.0f); - soundProcessor.MasterGain = 0.0f; //0.00001f; - soundProcessor.PressureGain = 0.1f; - soundProcessor.TurbulenceGain = 0.0f; - soundProcessor.Turbulence = 0.001f; - soundProcessor.SetAmbientPressure(ambientPressure); + // ---- Sound processor (stable version) ---- + soundProcessor = new SoundProcessor(sampleRate, pipeRadius * 2); + soundProcessor.Gain = 0.00001f; - Console.WriteLine("=== EngineScenario (torque‑driven RPM, throttle = pressure) ==="); - Console.WriteLine($"Crankshaft inertia: {crankshaft.Inertia}, friction: {crankshaft.FrictionConstant} + {crankshaft.FrictionViscous}*ω"); - Console.WriteLine($"Throttle range: {IdlePeakPressure/101325:F0} – {MaxPeakPressure/101325:F0} bar"); + // ---- Reverb ---- + reverb = new OutdoorExhaustReverb(sampleRate); + // Church: vast, highly reflective, bright + reverb.DryMix = 1.0f; // always full dry signal + reverb.EarlyMix = 0.5f; // distinct early reflections from distant walls + reverb.TailMix = 0.9f; // huge tail, almost as loud as the dry sound + reverb.Feedback = 0.9f; // long decay – roughly 3 s reverb time (with current delay lengths) + reverb.DampingFreq = 6000f; // bright: high‑frequency energy stays for a long time + reverb.MatrixCoeff = 0.5f; // default orthogonal mix + + Console.WriteLine("=== EngineScenario (Stable) ==="); + Console.WriteLine($"Crankshaft inertia: {crankshaft.Inertia}"); Console.WriteLine($"Pipe: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz"); } public override float Process() { - // 1. Map throttle to target peak pressure - double targetPressure = IdlePeakPressure + Throttle * (MaxPeakPressure - IdlePeakPressure); + // ---- RPM governor: adjust throttle to maintain idle when no user input ---- + double currentRPM = crankshaft.AngularVelocity * 60.0 / (2.0 * Math.PI); + double throttle = Math.Clamp(Throttle, 0.05, 1.0); // never let it drop below a tiny value + + // ---- Target combustion pressure ---- + double targetPressure = throttle * FullLoadPeakPressure; engineCyl.TargetPeakPressure = targetPressure; - // 2. Step the cylinder (adds torque to crankshaft, updates valve) + // ---- Simulate one timestep ---- engineCyl.Step(dt); - - // 3. Integrate crankshaft (applies friction, updates RPM) crankshaft.Step(dt); - - // 4. Set orifice area for coupling coupling.OrificeArea = engineCyl.OrificeArea; + solver.Step(); - // 5. Fluid solver step - float massFlow = solver.Step(); - float endPressure = (float)exhaustPipe.GetCellPressure(exhaustPipe.GetCellCount() - 1); + // ---- Update exit port with safety clamps ---- + UpdateExitPort(); - // 6. Audio - float audioSample = soundProcessor.Process(massFlow, endPressure); + // ---- Generate audio ---- + float dry = soundProcessor.Process(exitPort); + float wet = reverb.Process(dry); time += dt; stepCount++; - if (stepCount % LogInterval == 0) { - Console.WriteLine(audioSample); - } - if (stepCount % 1000 == 0 && false) - { - Console.WriteLine($"{time,5:F3} {crankshaft.AngularVelocity*60/(2*Math.PI),5:F0} RPM " + - $"Thr:{Throttle:F2} P_target:{targetPressure/101325:F1} bar " + - $"mflow:{massFlow,14:E4} Comb#{engineCyl.CombustionCount} Mis#{engineCyl.MisfireCount}"); - } - - return audioSample; + return wet; + } + + private void UpdateExitPort() + { + int last = exhaustPipe.GetCellCount() - 1; + double p = exhaustPipe.GetCellPressure(last); + double rho = exhaustPipe.GetCellDensity(last); + double vel = exhaustPipe.GetCellVelocity(last); + + // Clamp density to physically possible values + if (rho < 0.01) rho = 0.01; // never let it approach zero + if (rho > 50.0) rho = 50.0; // never let it become absurd + + // Clamp velocity to ± 500 m/s (safe subsonic) + vel = Math.Clamp(vel, -500.0, 500.0); + + double outflowMassFlow = rho * vel * pipeArea; + + // Clamp exit pressure to sensible range (0.1 – 20 bar) + p = Math.Clamp(p, 1e4, 2e6); + + exitPort.Pressure = p; + exitPort.Density = rho; + exitPort.Temperature = p / (rho * 287.05); + exitPort.MassFlowRate = -outflowMassFlow; + exitPort.SpecificEnthalpy = 0.0; } - // ---- Drawing (unchanged) ---- public override void Draw(RenderWindow target) { float winW = target.GetView().Size.X;