로켓 배기 플룸 쇼크 다이아몬드(Shock Diamond) & 1,000Hz 0-Alloc C# 마하 디스크 경사 충격파 해석 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Education.StemAero.ShockDiamond
{
    public readonly record struct RocketNozzleExitState(
        float AltitudeMeters,       // 비행 고도 (m)
        float NozzleExitDiameterM,  // 노즐 출구 직경 (m)
        float ExitMachNumber,       // 출구 마하수 (예: 3.2)
        float ExitPressureKpa,      // 출구 정압 (kPa, 예: 70 kPa)
        float ChamberTotalTempK     // 연소실 총온도 (K, 예: 3400 K)
    );

    public ref struct ShockDiamondTelemetry
    {
        public float AmbientPressureKpa;
        public float PressureRatioPePa;
        public float DiamondCellLengthMeters;
        public float EstimatedMachDiskDiameterMeters;
        public float PostShockTemperatureKelvin;
        public int VisibleDiamondCount;
    }

    public static class ShockDiamondEngineCore
    {
        private const float SeaLevelPressureKpa = 101.325f;
        private const float ScaleHeightMeters = 8400f; // 지구 대기 스케일 높이
        private const float SpecificHeatRatioGamma = 1.22f; // 로켓 연소가스 비열비

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void CalculatePlumeStructure(
            in RocketNozzleExitState state,
            ref ShockDiamondTelemetry telemetry)
        {
            // 1. 고도별 표준 대기압 계산: P_amb = P0 * exp(-h / H) (0-Alloc)
            float pAmbient = SeaLevelPressureKpa * MathF.Exp(-state.AltitudeMeters / ScaleHeightMeters);
            pAmbient = MathF.Max(pAmbient, 0.01f); // 진공 한계
            telemetry.AmbientPressureKpa = pAmbient;

            float pr = state.ExitPressureKpa / pAmbient;
            telemetry.PressureRatioPePa = pr;

            // 2. 쇼크 다이아몬드 셀 주기 길이 산출 L = d * sqrt(M^2 - 1) * sqrt(Pe/Pa)
            float machFactor = MathF.Sqrt(MathF.Max(0.1f, state.ExitMachNumber * state.ExitMachNumber - 1.0f));
            float expansionFactor = MathF.Sqrt(MathF.Max(0.1f, pr));
            float cellLength = state.NozzleExitDiameterM * machFactor * expansionFactor;
            telemetry.DiamondCellLengthMeters = cellLength;

            // 3. 마하 디스크 직경 및 충격파 후단 온도 추정
            if (pr < 0.7f) // 강한 과대팽창 (해수면 부근)
            {
                telemetry.EstimatedMachDiskDiameterMeters = state.NozzleExitDiameterM * (1.0f - pr) * 0.6f;
                // 수직 충격파 랭킨-위고니오 온도 상승
                float m1Sq = state.ExitMachNumber * state.ExitMachNumber;
                float tempRatio = (1f + 0.2f * m1Sq); // 정체온도 회복
                telemetry.PostShockTemperatureKelvin = state.ChamberTotalTempK * 0.85f;
                telemetry.VisibleDiamondCount = Math.Clamp((int)(4.5f / MathF.Max(0.5f, cellLength)), 1, 8);
            }
            else
            {
                telemetry.EstimatedMachDiskDiameterMeters = 0f; // 정규 반사 혹은 과소팽창
                telemetry.PostShockTemperatureKelvin = state.ChamberTotalTempK * 0.45f;
                telemetry.VisibleDiamondCount = (pr > 3.0f) ? 1 : 4; // 고고도에서는 단일 거대 플룸으로 팽창
            }
        }
    }
}