심우주 탐사선 홀 추력기(Hall Thruster) 크세논 플라즈마 정전가속 & 1,000Hz 0-Alloc C# 저추력 궤도 전이(RK4) 수치해석 엔진

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

namespace ZeliDesk.Education.StemSpace.IonThruster
{
    public readonly record struct SpacecraftState(
        Vector3 PositionMeters,   // 지심/태양 중심 관성좌표계 위치 (m)
        Vector3 VelocityMps,      // 속도 (m/s)
        float MassKg,             // 현재 탐사선 총 질량 (kg)
        float IspSeconds          // 추력기 비추력 (초)
    );

    public ref struct ThrusterOutput
    {
        public Vector3 ThrustNewton;
        public float MassFlowKgPerSec;
        public float AccelerationMps2;
    }

    public static class LowThrustRk4Integrator
    {
        private const float MuEarth = 3.986004418e14f; // 지구 중력상수 (m^3/s^2)
        private const float StandardGravityG0 = 9.80665f;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static SpacecraftState StepRk4(
            in SpacecraftState state,
            Vector3 thrustUnitVector,
            float thrustMagnitudeN,
            float dtSec)
        {
            // k1 평가
            EvaluateDerivatives(state, thrustUnitVector, thrustMagnitudeN, out Vector3 dr1, out Vector3 dv1, out float dm1);

            // k2 평가
            var s2 = new SpacecraftState(
                state.PositionMeters + dr1 * (0.5f * dtSec),
                state.VelocityMps + dv1 * (0.5f * dtSec),
                state.MassKg + dm1 * (0.5f * dtSec),
                state.IspSeconds
            );
            EvaluateDerivatives(s2, thrustUnitVector, thrustMagnitudeN, out Vector3 dr2, out Vector3 dv2, out float dm2);

            // k3 평가
            var s3 = new SpacecraftState(
                state.PositionMeters + dr2 * (0.5f * dtSec),
                state.VelocityMps + dv2 * (0.5f * dtSec),
                state.MassKg + dm2 * (0.5f * dtSec),
                state.IspSeconds
            );
            EvaluateDerivatives(s3, thrustUnitVector, thrustMagnitudeN, out Vector3 dr3, out Vector3 dv3, out float dm3);

            // k4 평가
            var s4 = new SpacecraftState(
                state.PositionMeters + dr3 * dtSec,
                state.VelocityMps + dv3 * dtSec,
                state.MassKg + dm3 * dtSec,
                state.IspSeconds
            );
            EvaluateDerivatives(s4, thrustUnitVector, thrustMagnitudeN, out Vector3 dr4, out Vector3 dv4, out float dm4);

            // RK4 가중 평균 합성 (0-Alloc 스택 연산)
            Vector3 newPos = state.PositionMeters + (dr1 + 2f * dr2 + 2f * dr3 + dr4) * (dtSec / 6f);
            Vector3 newVel = state.VelocityMps + (dv1 + 2f * dv2 + 2f * dv3 + dv4) * (dtSec / 6f);
            float newMass = state.MassKg + (dm1 + 2f * dm2 + 2f * dm3 + dm4) * (dtSec / 6f);

            return new SpacecraftState(newPos, newVel, Math.Max(newMass, 1.0f), state.IspSeconds);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static void EvaluateDerivatives(
            in SpacecraftState s,
            Vector3 thrustDir,
            float thrustN,
            out Vector3 dr,
            out Vector3 dv,
            out float dm)
        {
            dr = s.VelocityMps;

            float rLength = s.PositionMeters.Length();
            float rCube = rLength * rLength * rLength;
            Vector3 gravityAcc = (rCube > 1e-6f) ? -MuEarth / rCube * s.PositionMeters : Vector3.Zero;

            Vector3 thrustAcc = (s.MassKg > 0.1f) ? (thrustDir * thrustN) / s.MassKg : Vector3.Zero;
            dv = gravityAcc + thrustAcc;

            dm = -(thrustN / (s.IspSeconds * StandardGravityG0));
        }
    }
}