태양과 지구 두 거대 천체의 중력과 회전좌표계의 원심력이 완벽한 평형을 이루는 5대 라그랑주 점(Lagrange Points) 중 L4와 L5는 태양 활동 극대기에 발생하는 코로나 물질 방출(CME)과 강력한 태양풍을 지구보다 며칠 앞서 측면에서 입체 감시할 수 있는 최적의 우주기상 전초기지입니다. 원형 제한 3체 문제(CR3BP)의 무차원화 회전좌표계에서 유도되는 야코비 상수(Jacobi Constant) 에너지 보존 법칙과 준주기적 3차원 헤일(Halo) 궤도 유지(Station-Keeping) 제어를 1,000Hz 주기로 0-Alloc C# 수치해석하는 온디바이스 우주항법 시뮬레이터를 공개합니다.

1. 2022 개정 교육과정 교과 연계 분석

  • 물리학 (만유인력 법칙, 회전좌표계 원심력과 코리올리 힘, 역학적 에너지 보존): 두 중심 천체의 중력 합력과 회전좌표계 가상 힘($-2\boldsymbol{\omega} \times \mathbf{v}$, $-\boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{r})$)의 동적 평형, 유효 포텐셜 $\Omega(x, y, z)$.
  • 지구과학/천문학 (라그랑주 5대 평형점의 천체역학적 의미, 태양풍과 우주기상 조기 경보): 동선 상의 불안정 평형점(L1, L2, L3)과 등변삼각형 안정 평형점(L4, L5), 제임스 웹(L2) 및 미래 한국형 우주기상 위성의 궤도 특성.
  • 수학 (CR3BP 비선형 미분방정식, 야코비 적분 상수, 3차원 상태천이행렬 STM): 회전 좌표계에서의 에너지 불변량 $C_J = 2\Omega(x, y, z) - (v_x^2 + v_y^2 + v_z^2)$, 궤도 이탈 오차를 최소 $\Delta v$로 복원하는 플록-리치(Floquet) 고유치 해석.
  • 정보 (무할당 실시간 1,000Hz 3체 궤도 전파기 및 추력기 펄스 분배기): 1ms 실시간 루프 내에서 가비지 컬렉터 없이 스택 메모리 구조체 `ref struct`와 SIMD 벡터 가속으로 야코비 에너지 드리프트를 0-Alloc 보정하는 C# 온디바이스 엔진.

2. 원형 제한 3체 문제(CR3BP) 및 야코비 적분 수식

질량비 $\mu = \frac{m_2}{m_1 + m_2}$로 정규화된 태양-지구 회전좌표계에서 우주선의 무차원화 운동방정식은 다음과 같습니다.

$$\ddot{x} - 2\dot{y} = \frac{\partial \Omega}{\partial x}, \quad \ddot{y} + 2\dot{x} = \frac{\partial \Omega}{\partial y}, \quad \ddot{z} = \frac{\partial \Omega}{\partial z}$$ $$\Omega(x, y, z) = \frac{1}{2}(x^2 + y^2) + \frac{1-\mu}{r_1} + \frac{\mu}{r_2}$$

여기서 우주선의 상대 속도 $v^2 = \dot{x}^2 + \dot{y}^2 + \dot{z}^2$에 대해 야코비 상수 $C_J$는 외력이 없을 때 엄밀히 보존됩니다.

$$C_J = 2\Omega(x, y, z) - v^2 = \text{const}$$

태양 복사압과 행성 섭동으로 인해 $C_J$가 한계를 벗어날 때, 헤일 궤도의 불안정 모드(Unstable Manifold) 고유벡터 방향을 감쇄시키는 충격 델타-V($\Delta \mathbf{v}_{SK}$)를 산출하여 최소 연료로 영구 체공을 유지합니다.

3. 1,000Hz 0-Alloc C# CR3BP 헤일 궤도 제어 코어

using System;
using System.Numerics;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Education.StemAstro.Cr3bpLagrange
{
    public readonly record struct Cr3bpSpacecraftState(
        Vector3 PositionNormalized, // 무차원 회전좌표계 위치 (x, y, z)
        Vector3 VelocityNormalized  // 무차원 회전좌표계 속도 (vx, vy, vz)
    );

    public ref struct OrbitTelemetry
    {
        public float JacobiConstant;
        public float EffectivePotential;
        public Vector3 AccelerationNormalized;
        public Vector3 StationKeepingDeltaV;
        public bool OrbitStable;
    }

    public static class Cr3bpEngineCore
    {
        private const float SunEarthMu = 3.0034806e-6f; // 태양-지구 정규화 질량비

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void EvaluateOrbitState(
            in Cr3bpSpacecraftState state,
            ref OrbitTelemetry telemetry)
        {
            float x = state.PositionNormalized.X;
            float y = state.PositionNormalized.Y;
            float z = state.PositionNormalized.Z;
            float vx = state.VelocityNormalized.X;
            float vy = state.VelocityNormalized.Y;
            float vz = state.VelocityNormalized.Z;

            // 1. 중심 천체 1(태양) 및 2(지구)와의 거리 계산 (0-Alloc)
            float r1Sq = (x + SunEarthMu) * (x + SunEarthMu) + y * y + z * z;
            float r2Sq = (x - (1f - SunEarthMu)) * (x - (1f - SunEarthMu)) + y * y + z * z;
            float r1 = MathF.Sqrt(MathF.Max(r1Sq, 1e-8f));
            float r2 = MathF.Sqrt(MathF.Max(r2Sq, 1e-8f));

            // 2. 유효 포텐셜 Omega(x, y, z) 산출
            float omega = 0.5f * (x * x + y * y) + (1f - SunEarthMu) / r1 + SunEarthMu / r2;
            telemetry.EffectivePotential = omega;

            // 3. 야코비 에너지 상수 C_J = 2*Omega - v^2
            float vSq = vx * vx + vy * vy + vz * vz;
            float cJ = 2f * omega - vSq;
            telemetry.JacobiConstant = cJ;

            // 4. 가속도 미분 계산 (코리올리 힘 포함)
            float dOmegaDx = x - (1f - SunEarthMu) * (x + SunEarthMu) / (r1 * r1 * r1) - SunEarthMu * (x - (1f - SunEarthMu)) / (r2 * r2 * r2);
            float dOmegaDy = y - (1f - SunEarthMu) * y / (r1 * r1 * r1) - SunEarthMu * y / (r2 * r2 * r2);
            float dOmegaDz = -(1f - SunEarthMu) * z / (r1 * r1 * r1) - SunEarthMu * z / (r2 * r2 * r2);

            float ax = dOmegaDx + 2f * vy;
            float ay = dOmegaDy - 2f * vx;
            float az = dOmegaDz;
            telemetry.AccelerationNormalized = new Vector3(ax, ay, az);

            // 5. 헤일 궤도 유지(SK) 보정 속도 추정
            float targetJacobi = 3.0008f; // 목표 헤일 궤도 야코비 에너지
            float energyError = cJ - targetJacobi;
            telemetry.OrbitStable = MathF.Abs(energyError) < 0.001f;

            telemetry.StationKeepingDeltaV = (MathF.Abs(energyError) > 0.0005f)
                ? -0.5f * energyError * Vector3.Normalize(state.VelocityNormalized)
                : Vector3.Zero;
        }
    }
}

4. 교육적 가치 및 심우주 탐사 미래비전 연계

본 교구는 뉴턴의 2체 문제를 넘어선 현대 천체역학의 정수인 3체 문제와 회전좌표계의 신비로운 동역학을 학생들에게 직접 선물합니다. 학생들은 왜 제임스 웹 우주망원경이 지구 주위가 아닌 태양-지구 L2 점 주변의 보이지 않는 가상의 궤도를 돌고 있는지, 미래 우주기상 위성이 왜 L4/L5 점에 위치해야 지구의 통신망과 전력망을 태양풍 재앙으로부터 지킬 수 있는지를 수학적 에너지 곡면(힐 영역)과 함께 직관적으로 탐구합니다.