diff --git a/src/CounterDrone.Core/Algorithms/CloudExpansionModel.cs b/src/CounterDrone.Core/Algorithms/CloudExpansionModel.cs new file mode 100644 index 0000000..b35c34e --- /dev/null +++ b/src/CounterDrone.Core/Algorithms/CloudExpansionModel.cs @@ -0,0 +1,61 @@ +using System; +using CounterDrone.Core.Models; + +namespace CounterDrone.Core.Algorithms +{ + /// 云团膨胀模型——根据弹药和环境参数计算各阶段半径、密度、膨胀时间 + public class CloudExpansionModel + { + private readonly AmmunitionSpec _ammo; + private readonly CombatScene _env; + + public CloudExpansionModel(AmmunitionSpec ammo, CombatScene env) + { + _ammo = ammo; + _env = env; + } + + /// Phase 1 爆轰初始半径 (m) + public float InitialRadius => 3.3f * (float)Math.Pow(Math.Max(0.01, (float)_ammo.BurstChargeKg), 0.32f); + public float TurbulentExpansionK => (float)_ammo.TurbulentExpansionK; + + /// 指定时刻湍流膨胀半径 R(t) = R₀ + k√t + public float RadiusAt(float elapsedSeconds) + { + if (elapsedSeconds <= 0) return InitialRadius; + bool inPhase3 = elapsedSeconds > 30f; + if (!inPhase3) + return InitialRadius + TurbulentExpansionK * (float)Math.Sqrt(elapsedSeconds); + + // Phase 3: 高斯扩散 + float windSpeed = (float)_env.WindSpeed; + float x = Math.Max(1f, windSpeed * (elapsedSeconds - 30f)); + var cls = Kinematics.GetStabilityClass((WeatherType)_env.WeatherType, windSpeed); + float sY = Kinematics.SigmaY(cls, x); + float sZ = Kinematics.SigmaZ(cls, x); + float peakC = Kinematics.GaussianPeakConcentration((float)_ammo.SourceStrength, sY, sZ); + float effConc = (float)_ammo.EffectiveConcentration; + if (peakC <= effConc) return InitialRadius + TurbulentExpansionK * (float)Math.Sqrt(30f); + float sigma = (sY + sZ) / 2f; + float rPhase2 = InitialRadius + TurbulentExpansionK * (float)Math.Sqrt(30f); + return rPhase2 + sigma * (float)Math.Sqrt(2f * Math.Log(peakC / effConc)); + } + + /// 指定时刻中心浓度 + public float DensityAt(float elapsedSeconds) + { + float r = RadiusAt(elapsedSeconds); + float volume = (4f / 3f) * (float)Math.PI * r * r * r; + return volume > 0.001f ? (float)_ammo.SourceStrength / volume : (float)_ammo.CoreDensity; + } + + /// 达到指定半径所需时间 (s) + public float TimeToReach(float radius) + { + float r0 = InitialRadius; + if (radius <= r0) return 0; + float k = TurbulentExpansionK; + return (radius - r0) * (radius - r0) / (k * k); + } + } +} diff --git a/src/CounterDrone.Core/Algorithms/DefaultDefensePlanner.cs b/src/CounterDrone.Core/Algorithms/DefaultDefensePlanner.cs index 57d84a6..41d3861 100644 --- a/src/CounterDrone.Core/Algorithms/DefaultDefensePlanner.cs +++ b/src/CounterDrone.Core/Algorithms/DefaultDefensePlanner.cs @@ -349,41 +349,13 @@ namespace CounterDrone.Core.Algorithms private (float effectiveRadius, float expansionTime, float rPhase2) ComputeEffectiveRadius( DroneGroup threat, AmmunitionSpec ammo, CombatScene env) { - float totalFlightTime = threat.ArrivalTime * 2f; - float halfTime = totalFlightTime / 2f; - float windSpeed = (float)env.WindSpeed; - var weather = (WeatherType)env.WeatherType; - float avgSpeed = (float)threat.Target.TypicalSpeed / 3.6f; - - // Phase 1 初始半径(基于爆发药 TNT 当量) - float r0 = 3.3f * (float)Math.Pow(Math.Max(0.01, (float)ammo.BurstChargeKg), 0.32); - float k = (float)ammo.TurbulentExpansionK; - float maxDur = (float)ammo.MaxDuration; - - // Phase 2 湍流膨胀后的半径 - float phase2Time = Math.Min(halfTime, 30f); - float turbulentRadius = r0 + k * (float)Math.Sqrt(Math.Max(0, phase2Time)); - - float expansionTime = (float)Math.Pow((turbulentRadius - r0) / Math.Max(k, 0.01f), 2); - float effectiveR = turbulentRadius; - - // Phase 3 高斯扩散 - if (halfTime > 30f) - { - float x = Math.Max(1f, windSpeed * (halfTime - 30f)); - var cls = Kinematics.GetStabilityClass(weather, windSpeed); - float sY = Kinematics.SigmaY(cls, x); - float sZ = Kinematics.SigmaZ(cls, x); - float peakC = Kinematics.GaussianPeakConcentration((float)ammo.SourceStrength, sY, sZ); - float threshold = (float)ammo.EffectiveConcentration; - if (peakC > threshold) - { - float sigma = (sY + sZ) / 2f; - effectiveR = turbulentRadius + sigma * (float)Math.Sqrt(2f * Math.Log(peakC / threshold)); - } - } - - return (effectiveR, expansionTime, turbulentRadius); + var model = new CloudExpansionModel(ammo, env); + float tPhase2 = 30f; + float rPhase2 = model.RadiusAt(tPhase2); + float expansionTime = model.TimeToReach(rPhase2); + float halfTime = threat.ArrivalTime * 2f; // 总飞行时间的一半 + float effectiveR = model.RadiusAt(halfTime); + return (effectiveR, expansionTime, rPhase2); } private int CalcRoundsNeeded(DroneGroup threat, AmmunitionSpec ammo, @@ -399,14 +371,6 @@ namespace CounterDrone.Core.Algorithms return Math.Max(1, (int)Math.Ceiling(requiredCoverage / spacing)); - - var aerosolType = (AerosolType)ammo.AerosolType; - float neededExposure = aerosolType == AerosolType.ActiveMaterial ? 2f : 6f; - float singleCoverage = neededExposure * avgSpeed; - float requiredCoverage = singleCoverage; - - return Math.Max(1, - (int)Math.Ceiling((requiredCoverage - 2f * turbulentRadius) / spacing) + 1); } private float ComputeInterceptProbability(DroneGroup threat,