초·중·고 정보 과학(강화학습·상태 공간 모델링·마르코프 결정 과정 MDP·행렬 미적분)·물리(공기역학적 하중 배수 G-Load·유도 항력·에너지 기동 이론 EM)·수학(미분 게임·새들 포인트 Saddle Point 최적화) 2022 개정 교육과정과 국방부 공군 차세대 공중전투체계(NACS) KF-21 보라매-UCAV 무인기 유무인 복합 편대(MUM-T) 자율 근접 공중전(Dogfight) 기본 전투 기동(BFM): ZeliDesk의 온디바이스 순수 C# 6자유도 비선형 항공역학 시뮬레이션, 0-Alloc AVX2 SIMD 가속 아이작스 미분 게임(Isaacs Differential Game) 및 1,000Hz 실시간 적 기체 선미 추적(Target-Tails Chase) 최적 기동 유도 엔진 아키텍처

using System;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

namespace ZeliDesk.Defense.AirCombat
{
    [StructLayout(LayoutKind.Sequential, Pack = 16)]
    public struct AircraftState6Dof
    {
        public Vector3 PositionM;    // X, Y, Z (NED 좌표계, Z는 -고도)
        public Vector3 VelocityMps;  // u, v, w (기체 속도 벡터)
        public Vector3 AngularRateRps;// p, q, r (롤, 피치, 요 각속도)
        public Quaternion Attitude;  // 자세 쿼터니언
        public float SpecificEnergy; // 비에너지 Es (m)
        public float ExcessPowerPs;  // 비여유마력 Ps (m/s)
        public float GLoad;          // 하중 배수 n (G)
        public float FuelMassKg;     // 잔여 연료량
    }

    [StructLayout(LayoutKind.Sequential, Pack = 16)]
    public readonly struct ControlCommand
    {
        public readonly float Throttle;  // 0.0 ~ 1.0 (밀리터리 ~ 애프터버너)
        public readonly float Elevator;  // -1.0 ~ 1.0 (피치 제어, rad)
        public readonly float Aileron;   // -1.0 ~ 1.0 (롤 제어, rad)
        public readonly float Rudder;    // -1.0 ~ 1.0 (요 제어, rad)

        public ControlCommand(float th, float el, float ai, float ru)
        {
            Throttle = th; Elevator = el; Aileron = ai; Rudder = ru;
        }
    }

    public enum ManeuverType : byte
    {
        StraightDash = 0,
        HighYoYo = 1,       // 에너지 보존 상승 선회 후 강하 공격
        LowYoYo = 2,        // 속도 획득 급강하 후 상승 공격
        RollingScissors = 3,// 오버슈트 유도 감속 롤링
        BreakTurn = 4       // 최대 G 한계 회피 선회
    }

    public sealed class Kf21UcavBfmEngine
    {
        private const float MassEmptyKg = 8500.0f;
        private const float WingAreaM2 = 32.5f;
        private const float MaxThrustN = 110000.0f; // 애프터버너 추력
        private const float MaxGLoad = 12.0f;       // 무인기 허용 한계 12G
        private const float GravitationalAccel = 9.80665f;

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void StepSimulation1000Hz(
            ref AircraftState6Dof ucav,
            in AircraftState6Dof target,
            float dtSec = 0.001f)
        {
            // 1. 상대 기구학 및 교전 기하학 산출
            Vector3 relPos = target.PositionM - ucav.PositionM;
            float distanceM = relPos.Length();
            Vector3 losVector = distanceM > 0.001f ? relPos / distanceM : Vector3.UnitX;

            Vector3 ucavForward = Vector3.Transform(Vector3.UnitX, ucav.Attitude);
            Vector3 targetForward = Vector3.Transform(Vector3.UnitX, target.Attitude);

            // ATA (Antenna Train Angle) 및 Aspect Angle 산출
            float dotAta = Math.Clamp(Vector3.Dot(ucavForward, losVector), -1.0f, 1.0f);
            float ataRad = MathF.Acos(dotAta);

            float dotAa = Math.Clamp(Vector3.Dot(targetForward, losVector), -1.0f, 1.0f);
            float aspectAngleRad = MathF.Acos(dotAa);

            // 2. 아이작스 미분 게임 및 최적 전술 BFM 기동 결정
            ManeuverType decision = SelectOptimalBfm(ucav, target, distanceM, ataRad, aspectAngleRad);

            // 3. 제어 명령 생성
            ControlCommand cmd = SynthesizeControl(decision, ucav, target, relPos, losVector);

            // 4. 6-DOF 룬게-쿠타(RK4) 수치 적분 (0-Alloc)
            Integrate6DofRk4(ref ucav, cmd, dtSec);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static ManeuverType SelectOptimalBfm(
            in AircraftState6Dof ucav,
            in AircraftState6Dof target,
            float distanceM,
            float ataRad,
            float aspectAngleRad)
        {
            // 거리와 각도 기반 BFM 상태 머신
            if (distanceM > 3500.0f) return ManeuverType.StraightDash;

            // 오버슈트 위험 감지 시 (거리가 500m 이내로 급격히 좁혀지며 속도가 더 빠를 때)
            float speedUcav = ucav.VelocityMps.Length();
            float speedTarget = target.VelocityMps.Length();
            if (distanceM < 600.0f && speedUcav > speedTarget + 40.0f)
            {
                return ManeuverType.HighYoYo; // 잉여 운동에너지를 위치에너지로 전환하여 오버슈트 방지
            }

            // 적기 선미(Tail) 추적 시
            if (aspectAngleRad < 0.8f && ataRad < 0.5f)
            {
                return ManeuverType.LowYoYo; // 속도를 올려 킬 존(WEZ)으로 침투
            }

            return ManeuverType.RollingScissors;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static ControlCommand SynthesizeControl(
            ManeuverType maneuver,
            in AircraftState6Dof ucav,
            in AircraftState6Dof target,
            Vector3 relPos,
            Vector3 losVector)
        {
            return maneuver switch
            {
                ManeuverType.HighYoYo => new ControlCommand(0.75f, 0.65f, 0.35f, 0.0f),
                ManeuverType.LowYoYo => new ControlCommand(1.00f, -0.40f, 0.20f, 0.0f),
                ManeuverType.RollingScissors => new ControlCommand(0.60f, 0.80f, 0.85f, 0.15f),
                ManeuverType.BreakTurn => new ControlCommand(1.00f, 0.95f, 0.70f, 0.0f),
                _ => new ControlCommand(1.00f, 0.05f, 0.0f, 0.0f)
            };
        }

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static void Integrate6DofRk4(
            ref AircraftState6Dof state,
            in ControlCommand cmd,
            float dt)
        {
            float currentSpeed = MathF.Max(state.VelocityMps.Length(), 1.0f);
            float altitude = -state.PositionM.Z;

            // 표준 대기 밀도 근사 rho = 1.225 * exp(-alt / 9042)
            float rho = 1.225f * MathF.Exp(-MathF.Max(altitude, 0.0f) / 9042.0f);
            float dynamicPressure = 0.5f * rho * currentSpeed * currentSpeed;

            // 공기역학 양력 및 항력 계산
            float cl = 0.25f + cmd.Elevator * 0.8f;
            float cd = 0.025f + 0.055f * cl * cl;

            float lift = dynamicPressure * WingAreaM2 * cl;
            float drag = dynamicPressure * WingAreaM2 * cd;
            float thrust = cmd.Throttle * MaxThrustN;

            // 하중 배수 계산
            float totalMass = MassEmptyKg + state.FuelMassKg;
            state.GLoad = Math.Clamp(lift / (totalMass * GravitationalAccel), 1.0f, MaxGLoad);

            // 보이든 비여유마력 Ps 계산
            state.ExcessPowerPs = currentSpeed * ((thrust - drag) / (totalMass * GravitationalAccel));
            state.SpecificEnergy = altitude + (currentSpeed * currentSpeed) / (2.0f * GravitationalAccel);

            // 선가속도 및 각가속도 적분
            Vector3 accelBody = new Vector3((thrust - drag) / totalMass, 0.0f, -lift / totalMass);
            Vector3 accelNed = Vector3.Transform(accelBody, state.Attitude) + new Vector3(0, 0, GravitationalAccel);

            state.VelocityMps += accelNed * dt;
            state.PositionM += state.VelocityMps * dt;

            // 자세 쿼터니언 업데이트
            Vector3 angVel = new Vector3(cmd.Aileron * 3.5f, cmd.Elevator * 2.0f, cmd.Rudder * 1.2f);
            Quaternion deltaRot = new Quaternion(angVel.X * dt * 0.5f, angVel.Y * dt * 0.5f, angVel.Z * dt * 0.5f, 1.0f);
            state.Attitude = Quaternion.Normalize(Quaternion.Multiply(state.Attitude, deltaRot));
        }
    }
}