우주 이온추진기(Ion Thruster) 차일드-랭뮤어(Child-Langmuir) 공간전하 한계 & 1,000Hz 0-Alloc C# 가속그리드 침식 제어 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.ElectricPropulsion;

public readonly ref struct IonThrusterState
{
    public readonly double BeamCurrentAmperes;
    public readonly double ThrustMilliNewtons;
    public readonly double SpecificImpulseSec;
    public readonly double GridErosionRateNmPerHour;
    public readonly bool IsDirectImpingementRisk;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public IonThrusterState(double iBeam, double thrustMn, double isp, double erosionNmH, bool risk)
    {
        BeamCurrentAmperes = iBeam;
        ThrustMilliNewtons = thrustMn;
        SpecificImpulseSec = isp;
        GridErosionRateNmPerHour = erosionNmH;
        IsDirectImpingementRisk = risk;
    }
}

public ref struct IonThrusterGridSimulator
{
    private const double Epsilon0 = 8.8541878e-12;
    private const double ElectronCharge = 1.602176634e-19;
    private const double XenonMassKg = 2.18e-25; // Xe 원자 질량
    private const double GridGapMeters = 0.001; // 1mm 그리드 간극
    private const double GridAreaM2 = 0.07; // 30cm 직경 그리드 유효 면적
    private const double GravityG0 = 9.80665;

    private double _screenVoltageV;
    private double _accelVoltageV;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public IonThrusterGridSimulator(double screenV = 1100.0, double accelV = -250.0)
    {
        _screenVoltageV = screenV;
        _accelVoltageV = accelV;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public IonThrusterState Step(double requestedBeamCurrentA, double neutralGasLeakFraction = 0.15)
    {
        double vTotal = _screenVoltageV + Math.Abs(_accelVoltageV);

        // 1. 차일드-랭뮤어 공간전하 제한 전류 밀도 산출 (A/m^2)
        double prefactor = (4.0 * Epsilon0 / 9.0) * Math.Sqrt(2.0 * ElectronCharge / XenonMassKg);
        double jChildLangmuir = prefactor * Math.Pow(vTotal, 1.5) / (GridGapMeters * GridGapMeters);
        double maxBeamCurrentA = jChildLangmuir * GridAreaM2 * 0.7; // 그리드 개구율 70%

        // 2. 유효 빔 전류 및 펄비언스 불일치(Direct Impingement) 위험 판정
        double actualBeamCurrent = Math.Min(requestedBeamCurrentA, maxBeamCurrentA);
        double perveanceRatio = actualBeamCurrent / Math.Max(0.01, maxBeamCurrentA);
        bool directImpingement = perveanceRatio > 0.88 || perveanceRatio < 0.35; // 과소/과대 시 빔 발산 충돌

        // 3. 이온 배출 속도, 비추력 및 추력 산출
        double ionVelocityMps = Math.Sqrt(2.0 * ElectronCharge * _screenVoltageV / XenonMassKg);
        double specificImpulseSec = (ionVelocityMps * 0.98) / GravityG0; // 발산 손실 보정
        double massFlowRateKgS = (actualBeamCurrent / ElectronCharge) * XenonMassKg * (1.0 + neutralGasLeakFraction);
        double thrustNewtons = (actualBeamCurrent / ElectronCharge) * XenonMassKg * ionVelocityMps * 0.98;
        double thrustMilliNewtons = thrustNewtons * 1e3;

        // 4. 전하교환(CEX) 이온에 의한 몰리브덴 가속 그리드 스퍼터 마모율 (nm/hour)
        double accelEnergyEv = Math.Abs(_accelVoltageV);
        double sputterYieldMo = 0.012 * Math.Sqrt(Math.Max(0.0, (accelEnergyEv - 40.0) / 100.0)); // 40eV 문턱 전위
        double cexCurrentA = actualBeamCurrent * neutralGasLeakFraction * 0.08;
        double erosionRateMps = (cexCurrentA * sputterYieldMo * 0.096) / (10280.0 * ElectronCharge * 6.022e23); // Mo 밀도 10280 kg/m^3
        double erosionRateNmPerHour = erosionRateMps * 1e9 * 3600.0;

        return new IonThrusterState(actualBeamCurrent, thrustMilliNewtons, specificImpulseSec, erosionRateNmPerHour, directImpingement);
    }
}