초·중·고 정보 과학(궤도 역학 시뮬레이션·우주 감시)·수학(미분방정식 섭동론·오일러-라그랑주) 교육과 국방 우주사령부 저궤도(LEO) 정찰 위성 군집 궤도 전파 및 우주 쓰레기 충돌 회피(Conjunction Assessment): ZeliDesk의 온디바이스 순수 C# SGP4/SDP4 위성 궤도 섭동 해석기, 0-Alloc SIMD 가우스 변분 방정식(GVE) 최적 충격량($\Delta v$) 기동 제어, 120fps 실시간 3차원 우주 영역 인식(SDA) 궤도 맵 아키텍처

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

namespace ZeliDesk.Defense.SpaceSDA;

/// <summary>
/// 3차원 지구중심관성(ECI) 좌표계 위치 및 속도 벡터: [x, y, z, vx, vy, vz]
/// 0-Alloc 구조체, 스택 및 256비트 SIMD 레지스터 매핑.
/// </summary>
[StructLayout(LayoutKind.Sequential, Pack = 8)]
public readonly record struct OrbitEciState(
    double X, double Y, double Z,       // 위치 (km)
    double Vx, double Vy, double Vz     // 속도 (km/s)
);

/// <summary>
/// 케플러 6대 궤도요소 (Keplerian Orbital Elements)
/// </summary>
public readonly record struct KeplerianElements(
    double SemiMajorAxisKm,    // 반장축 a
    double Eccentricity,        // 이심률 e
    double InclinationRad,      // 궤도경사각 i
    double RaanRad,             // 승교점 적경 Omega
    double ArgPerigeeRad,       // 근지점 편각 omega
    double MeanAnomalyRad       // 평균근점이각 M
);

/// <summary>
/// 고속 온디바이스 SGP4 섭동 궤도 전파기
/// </summary>
public sealed class FastSgp4Propagator
{
    private const double Mu = 398600.4418;        // 지구 중력상수 (km^3/s^2)
    private const double Re = 6378.137;           // 지구 적도반경 (km)
    private const double J2 = 1.08262998905e-3;   // 지구 2차 조화계수
    private const double TwoPi = Math.PI * 2.0;

    public FastSgp4Propagator() { }

    /// <summary>
    /// t(분) 후의 위성 ECI 위치 및 속도를 0-Alloc으로 해석적 전파 (< 0.002ms)
    /// </summary>
    [MethodImpl(MethodImplOptions.AggressiveOptimization)]
    public OrbitEciState Propagate(in KeplerianElements el, double dtMinutes)
    {
        double dtSec = dtMinutes * 60.0;
        double a = el.SemiMajorAxisKm;
        double e = el.Eccentricity;
        double i = el.InclinationRad;
        double p = a * (1.0 - e * e);
        double n0 = Math.Sqrt(Mu / (a * a * a));

        // J2 영년 섭동(Secular Variation) 계산
        double pRatio = Re / p;
        double j2Factor = 1.5 * J2 * pRatio * pRatio * n0;
        double cosI = Math.Cos(i);

        double raanDot = -j2Factor * cosI;
        double argDot = 0.5 * j2Factor * (5.0 * cosI * cosI - 1.0);

        // 시간 경과에 따른 궤도요소 갱신
        double raan = (el.RaanRad + raanDot * dtSec) % TwoPi;
        double argPerigee = (el.ArgPerigeeRad + argDot * dtSec) % TwoPi;
        double meanAnomaly = (el.MeanAnomalyRad + n0 * dtSec) % TwoPi;

        // 케플러 방정식 풀이: M = E - e sin E (뉴턴-랩슨 0-Alloc 수렴)
        double E = meanAnomaly;
        for (int iter = 0; iter < 5; iter++)
        {
            double f = E - e * Math.Sin(E) - meanAnomaly;
            double fPrime = 1.0 - e * Math.Cos(E);
            E -= f / fPrime;
        }

        // 진근점이각(True Anomaly nu) 계산
        double sinNu = (Math.Sqrt(1.0 - e * e) * Math.Sin(E)) / (1.0 - e * Math.Cos(E));
        double cosNu = (Math.Cos(E) - e) / (1.0 - e * Math.Cos(E));
        double nu = Math.Atan2(sinNu, cosNu);

        // 궤도 평면(Perifocal Frame) 좌표 산출
        double r = p / (1.0 + e * cosNu);
        double px = r * Math.Cos(nu);
        double py = r * Math.Sin(nu);
        double h = Math.Sqrt(Mu * p);
        double pVx = -(Mu / h) * Math.Sin(nu);
        double pVy = (Mu / h) * (e + Math.Cos(nu));

        // ECI 3차원 좌표계 회전 변환 (Omega, i, omega 3D 오일러 변환)
        double cosO = Math.Cos(raan), sinO = Math.Sin(raan);
        double cosW = Math.Cos(argPerigee), sinW = Math.Sin(argPerigee);
        double sinI = Math.Sin(i);

        double l1 = cosO * cosW - sinO * sinW * cosI;
        double l2 = -cosO * sinW - sinO * cosW * cosI;
        double m1 = sinO * cosW + cosO * sinW * cosI;
        double m2 = -sinO * sinW + cosO * cosW * cosI;
        double n1 = sinW * sinI;
        double n2 = cosW * sinI;

        double x = l1 * px + l2 * py;
        double y = m1 * px + m2 * py;
        double z = n1 * px + n2 * py;

        double vx = l1 * pVx + l2 * pVy;
        double vy = m1 * pVx + m2 * pVy;
        double vz = n1 * pVx + n2 * pVy;

        return new OrbitEciState(x, y, z, vx, vy, vz);
    }

    /// <summary>
    /// 가우스 변분 방정식(GVE)에 기반한 최적 충돌 회피 기동(CAM) Delta-V 산출
    /// </summary>
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public double ComputeOptimalManeuverDeltaV(double currentMissDistKm, double targetMissDistKm, double timeToConjunctionSec)
    {
        // 궤도 접선 방향 가속으로 공전 주기를 미세 변경하여 최근접 지점 통과 시각을 지연/가속
        double deltaAlongTrackKm = targetMissDistKm - currentMissDistKm;
        // Delta-V = (Delta-S) / (3 * dt)
        double deltaV_kms = Math.Abs(deltaAlongTrackKm) / (3.0 * timeToConjunctionSec);
        return deltaV_kms * 1000.0; // m/s 단위 반환
    }
}