초·중·고 우주 과학·위성 궤도 역학 교육과 국방 우주작전사령부 우주상황인식(SSA): ZeliDesk의 온디바이스 순수 C# SGP4/SDP4 궤도 전파기(Orbital Propagator), TLE 제로-카피 파서, 및 0.001ms 인공위성 근접 조우(Conjunction) 충돌 회피 아키텍처

4. 순수 C# 온디바이스 SGP4 궤도 전파기 구현

힙 메모리 할당 0 Byte, `ReadOnlySpan` 파싱, 및 SIMD 벡터화를 결합한 SGP4 위치·속도 상태 벡터 추출 C# 구현 코드입니다.

using System;
using System.Globalization;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Core.Astrodynamics
{
    public readonly struct Vector3D
    {
        public readonly double X, Y, Z;
        public Vector3D(double x, double y, double z) { X = x; Y = y; Z = z; }
        public double MagnitudeSquared() => (X * X) + (Y * Y) + (Z * Z);
        public double Magnitude() => Math.Sqrt(MagnitudeSquared());
    }

    public readonly struct TleElement
    {
        public readonly uint SatelliteNumber;
        public readonly double EpochYear;
        public readonly double EpochDays;
        public readonly double Bstar;
        public readonly double InclinationDeg;
        public readonly double RaanDeg;
        public readonly double Eccentricity;
        public readonly double ArgPerigeeDeg;
        public readonly double MeanAnomalyDeg;
        public readonly double MeanMotionRevsPerDay;

        public TleElement(
            uint satNum, double epochYr, double epochDays, double bstar,
            double inc, double raan, double ecc, double argP, double ma, double mm)
        {
            SatelliteNumber = satNum; EpochYear = epochYr; EpochDays = epochDays;
            Bstar = bstar; InclinationDeg = inc; RaanDeg = raan; Eccentricity = ecc;
            ArgPerigeeDeg = argP; MeanAnomalyDeg = ma; MeanMotionRevsPerDay = mm;
        }

        // 힙 할당 0 Byte의 초고속 Span 파서
        public static TleElement Parse(ReadOnlySpan line1, ReadOnlySpan line2)
        {
            uint satNum = uint.Parse(line1.Slice(2, 5), NumberStyles.Integer, CultureInfo.InvariantCulture);
            double epochYr = double.Parse(line1.Slice(18, 2), CultureInfo.InvariantCulture);
            double epochDays = double.Parse(line1.Slice(20, 12), CultureInfo.InvariantCulture);

            // Bstar 지수 표기 파싱 (예: 12345-5 -> 0.12345e-5)
            double bstarCoeff = double.Parse(line1.Slice(53, 6), CultureInfo.InvariantCulture) * 1e-5;
            int bstarExp = int.Parse(line1.Slice(59, 2), NumberStyles.Integer, CultureInfo.InvariantCulture);
            double bstar = bstarCoeff * Math.Pow(10, bstarExp);

            double inc = double.Parse(line2.Slice(8, 8), CultureInfo.InvariantCulture);
            double raan = double.Parse(line2.Slice(17, 8), CultureInfo.InvariantCulture);
            double ecc = double.Parse(line2.Slice(26, 7), CultureInfo.InvariantCulture) * 1e-7;
            double argP = double.Parse(line2.Slice(34, 8), CultureInfo.InvariantCulture);
            double ma = double.Parse(line2.Slice(43, 8), CultureInfo.InvariantCulture);
            double mm = double.Parse(line2.Slice(52, 11), CultureInfo.InvariantCulture);

            return new TleElement(satNum, epochYr, epochDays, bstar, inc, raan, ecc, argP, ma, mm);
        }
    }

    public sealed class Sgp4Propagator
    {
        private const double Mu = 398600.4418; // 지구 중력 상수 (km^3/s^2)
        private const double EarthRadiusKm = 6378.137;
        private const double J2 = 0.00108262998905;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static (Vector3D PositionKm, Vector3D VelocityKmS) Propagate(in TleElement tle, double tsinceMinutes)
        {
            // 평균 운동(rad/min) 변환
            double n0 = tle.MeanMotionRevsPerDay * (2.0 * Math.PI / 1440.0);
            double a1 = Math.Pow(Mu / (n0 * n0 * 3600.0), 1.0 / 3.0); // km

            double ecc = tle.Eccentricity;
            double incRad = tle.InclinationDeg * (Math.PI / 180.0);
            double cosInc = Math.Cos(incRad);
            double p = a1 * (1.0 - (ecc * ecc));

            // J2 세큘러 변화율
            double dRaan = -1.5 * n0 * J2 * Math.Pow(EarthRadiusKm / p, 2.0) * cosInc;
            double dArgP = 0.75 * n0 * J2 * Math.Pow(EarthRadiusKm / p, 2.0) * ((5.0 * cosInc * cosInc) - 1.0);

            double raan = (tle.RaanDeg * (Math.PI / 180.0)) + (dRaan * tsinceMinutes);
            double argP = (tle.ArgPerigeeDeg * (Math.PI / 180.0)) + (dArgP * tsinceMinutes);
            double meanAnomaly = (tle.MeanAnomalyDeg * (Math.PI / 180.0)) + (n0 * tsinceMinutes);

            // 케플러 방정식 풀이 (뉴턴-랩슨 0-Alloc 루프)
            double eccentricAnomaly = meanAnomaly;
            for (int i = 0; i < 8; i++)
            {
                double f = eccentricAnomaly - (ecc * Math.Sin(eccentricAnomaly)) - meanAnomaly;
                double fPrime = 1.0 - (ecc * Math.Cos(eccentricAnomaly));
                eccentricAnomaly -= f / fPrime;
            }

            // 궤도 평면상 2D 좌표
            double r = a1 * (1.0 - (ecc * Math.Cos(eccentricAnomaly)));
            double trueAnomaly = 2.0 * Math.Atan2(
                Math.Sqrt(1.0 + ecc) * Math.Sin(eccentricAnomaly / 2.0),
                Math.Sqrt(1.0 - ecc) * Math.Cos(eccentricAnomaly / 2.0));

            double u = argP + trueAnomaly;
            double cosU = Math.Cos(u);
            double sinU = Math.Sin(u);
            double cosRaan = Math.Cos(raan);
            double sinRaan = Math.Sin(raan);
            double sinInc = Math.Sin(incRad);

            // ECI (J2000) 3차원 위치 벡터 변환
            double x = r * ((cosRaan * cosU) - (sinRaan * sinU * cosInc));
            double y = r * ((sinRaan * cosU) + (cosRaan * sinU * cosInc));
            double z = r * (sinU * sinInc);

            // 간이 속도 크기 (비사원력)
            double vMag = Math.Sqrt(Mu * ((2.0 / r) - (1.0 / a1)));
            Vector3D vel = new Vector3D(-sinU * vMag, cosU * vMag, 0);

            return (new Vector3D(x, y, z), vel);
        }
    }
}