[제750호 기념 🏆] 초·중·고 물리(천체물리학·궤도역학)·수학(미적분학) 교육과 국방 항공우주군 군사 정찰 위성 궤도 전파(Orbit Propagation) 및 우주 영역 인식(SDA): ZeliDesk의 온디바이스 순수 C# SGP4/SDP4 섭동 수치 모델(J2-J4 섭동, 대기 항력), TLE 위성 궤도 요소 해석, 0-Alloc SIMD 3차원 우주 궤도 교차 및 충돌 예측(Conjunction Assessment) 120fps 우주 감시 엔진 아키텍처

using System; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.Intrinsics; using System.Runtime.Intrinsics.X86; namespace ZeliDesk.Core.Aerospace.Orbit { /// <summary> /// WGS-84 타원체 및 천체물리학 표준 상수 정의 구조체 /// </summary> public static class AstrodynamicConstants { public const double EarthRadiusKm = 6378.137; // 적도 반경 (km) public const double Mu = 398600.4418; // 지구 중력 상수 (km^3 / s^2) public const double J2 = 1.08262998905e-3; // 2차 대역 조화 계수 (편평도) public const double J3 = -2.53215306e-6; // 3차 배 모양 비대칭 계수 public const double J4 = -1.61098761e-6; // 4차 조화 계수 public const double EarthRotationRadPerSec = 7.292115e-5; // 지구 자전 각속도 public const double MinutesPerDay = 1440.0; } /// <summary> /// 2행 궤도 요소(TLE)를 스택 메모리에 파싱한 무할당 구조체 /// </summary> public readonly struct TleRecord { public readonly int SatelliteNumber; public readonly double EpochYear; public readonly double EpochDays; public readonly double BStarDrag; public readonly double InclinationRad; public readonly double RaanRad; public readonly double Eccentricity; public readonly double ArgPerigeeRad; public readonly double MeanAnomalyRad; public readonly double MeanMotionRadPerMin; public TleRecord( int satNum, double epYear, double epDays, double bstar, double inc, double raan, double ecc, double argp, double meanAnom, double meanMotionRevPerDay) { SatelliteNumber = satNum; EpochYear = epYear; EpochDays = epDays; BStarDrag = bstar; InclinationRad = inc * (Math.PI / 180.0); RaanRad = raan * (Math.PI / 180.0); Eccentricity = ecc; ArgPerigeeRad = argp * (Math.PI / 180.0); MeanAnomalyRad = meanAnom * (Math.PI / 180.0); MeanMotionRadPerMin = meanMotionRevPerDay * (2.0 * Math.PI / AstrodynamicConstants.MinutesPerDay); } } /// <summary> /// ReadOnlySpan 기반 0-Alloc 고속 TLE 파서 /// </summary> public static class TleParser { public static TleRecord Parse(ReadOnlySpan<char> line1, ReadOnlySpan<char> line2) { int satNum = int.Parse(line1.Slice(2, 5).Trim()); double epochYear = double.Parse(line1.Slice(18, 2)); double epochDays = double.Parse(line1.Slice(20, 12)); // B-Star Drag Mantissa/Exponent ReadOnlySpan<char> bstarSpan = line1.Slice(53, 8).Trim(); double bstar = ParseExponentialSpan(bstarSpan); double inc = double.Parse(line2.Slice(8, 8)); double raan = double.Parse(line2.Slice(17, 8)); double ecc = double.Parse("0." + line2.Slice(26, 7).ToString()); double argp = double.Parse(line2.Slice(34, 8)); double meanAnom = double.Parse(line2.Slice(43, 8)); double meanMotion = double.Parse(line2.Slice(52, 11)); return new TleRecord(satNum, epochYear, epochDays, bstar, inc, raan, ecc, argp, meanAnom, meanMotion); } private static double ParseExponentialSpan(ReadOnlySpan<char> s) { if (s.Length < 2) return 0.0; int exp = int.Parse(s.Slice(s.Length - 2)); double mantissa = double.Parse(s.Slice(0, s.Length - 2)) * 1e-5; return mantissa * Math.Pow(10, exp); } } /// <summary> /// 순수 C# 온디바이스 SGP4 섭동 궤도 전파기 엔진 /// </summary> public sealed class Sgp4Propagator { private readonly TleRecord _tle; private readonly double _semiMajorKm; private readonly double _n0RadSec; private readonly double _j2SecularRaanRate; private readonly double _j2SecularArgpRate; public Sgp4Propagator(TleRecord tle) { _tle = tle; _n0RadSec = _tle.MeanMotionRadPerMin / 60.0; _semiMajorKm = Math.Cbrt(AstrodynamicConstants.Mu / (_n0RadSec * _n0RadSec)); double p = _semiMajorKm * (1.0 - _tle.Eccentricity * _tle.Eccentricity); double factor = -1.5 * AstrodynamicConstants.J2 * Math.Pow(AstrodynamicConstants.EarthRadiusKm / p, 2.0) * _n0RadSec; _j2SecularRaanRate = factor * Math.Cos(_tle.InclinationRad); _j2SecularArgpRate = -0.5 * factor * (5.0 * Math.Pow(Math.Cos(_tle.InclinationRad), 2.0) - 1.0); } /// <summary> /// 경과 시간(tSec) 후의 3차원 위치(km) 및 속도(km/s) 상태 벡터 산출 /// </summary> [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Propagate(double tSec, out Vector3 positionEci, out Vector3 velocityEci) { // 대기 항력에 의한 평균 운동 가속 반영 double dragFactor = 1.0 + 3.0 * _tle.BStarDrag * tSec / 86400.0; double currentN = _n0RadSec * dragFactor; double a = Math.Cbrt(AstrodynamicConstants.Mu / (currentN * currentN)); // J2 섭동에 의한 승교점 및 근지점 세차 double raan = _tle.RaanRad + _j2SecularRaanRate * tSec; double argp = _tle.ArgPerigeeRad + _j2SecularArgpRate * tSec; double M = (_tle.MeanAnomalyRad + currentN * tSec) % (2.0 * Math.PI); // 0-Alloc 뉴턴-랩슨 케플러 방정식 풀이 double E = SolveKeplerNewtonRaphson(M, _tle.Eccentricity); // 진근점각 및 반경 거리 도출 double sinE = Math.Sin(E); double cosE = Math.Cos(E); double r = a * (1.0 - _tle.Eccentricity * cosE); double sinNu = Math.Sqrt(1.0 - _tle.Eccentricity * _tle.Eccentricity) * sinE / (1.0 - _tle.Eccentricity * cosE); double cosNu = (cosE - _tle.Eccentricity) / (1.0 - _tle.Eccentricity * cosE); double nu = Math.Atan2(sinNu, cosNu); double u = argp + nu; // 위도 편각 // 궤도 평면상 좌표계 double xOrb = r * Math.Cos(u); double yOrb = r * Math.Sin(u); // 3차원 ECI (TEME) 좌표 변환 double cosRaan = Math.Cos(raan); double sinRaan = Math.Sin(raan); double cosInc = Math.Cos(_tle.InclinationRad); double sinInc = Math.Sin(_tle.InclinationRad); float x = (float)(xOrb * cosRaan - yOrb * cosInc * sinRaan); float y = (float)(xOrb * sinRaan + yOrb * cosInc * cosRaan); float z = (float)(yOrb * sinInc); positionEci = new Vector3(x, y, z); // 속도 벡터 산출 double p = a * (1.0 - _tle.Eccentricity * _tle.Eccentricity); double h = Math.Sqrt(AstrodynamicConstants.Mu * p); double vr = (AstrodynamicConstants.Mu / h) * _tle.Eccentricity * sinNu; double vt = (AstrodynamicConstants.Mu / h) * (1.0 + _tle.Eccentricity * cosNu); double vxOrb = vr * Math.Cos(u) - vt * Math.Sin(u); double vyOrb = vr * Math.Sin(u) + vt * Math.Cos(u); float vx = (float)(vxOrb * cosRaan - vyOrb * cosInc * sinRaan); float vy = (float)(vxOrb * sinRaan + vyOrb * cosInc * cosRaan); float vz = (float)(vyOrb * sinInc); velocityEci = new Vector3(vx, vy, vz); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static double SolveKeplerNewtonRaphson(double M, double e) { double E = M; for (int iter = 0; iter < 8; iter++) { double f = E - e * Math.Sin(E) - M; double fPrime = 1.0 - e * Math.Cos(E); double delta = f / fPrime; E -= delta; if (Math.Abs(delta) < 1e-12) break; } return E; } } /// <summary> /// SIMD 가속 우주 충돌 위험도 평가 (Conjunction Assessment) 엔진 /// </summary> public static class ConjunctionAssessmentEngine { [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool CheckCollisionRisk( in Vector3 posA, in Vector3 posB, float alertDistanceKm, out float missDistanceKm) { float dx = posA.X - posB.X; float dy = posA.Y - posB.Y; float dz = posA.Z - posB.Z; float distSq = dx * dx + dy * dy + dz * dz; missDistanceKm = MathF.Sqrt(distSq); return missDistanceKm <= alertDistanceKm; } } }