초·중·고 물리·미적분학 교육과 국방 전술 포병 사격 통제(Fire Control System) 탄도학: ZeliDesk의 온디바이스 순수 C# 4차 룬게-쿠타(RK4) 수치 적분기, G1/G7 마하 항력 함수, 코리올리 힘(Coriolis Effect) 3차원 지구 자전 보정, 및 0.001mil 정밀도 실시간 탄착점 산출 아키텍처

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

namespace ZeliDesk.Defense.Ballistics
{
    public readonly struct BallisticState
    {
        public readonly Vector3 Position; // x: 사거리(m), y: 편각(m), z: 고도(m)
        public readonly Vector3 Velocity; // vx, vy, vz (m/s)

        public BallisticState(Vector3 pos, Vector3 vel)
        {
            Position = pos;
            Velocity = vel;
        }
    }

    /// 

    /// ZeliDesk 온디바이스 순수 C# 4차 룬게-쿠타(RK4) 탄도학 솔버
    /// 
    public sealed class BallisticSolverCore
    {
        private const float StandardGravity = 9.80665f;
        private const float EarthAngularSpeed = 7.292115e-5f; // rad/s

        /// 

        /// 대기 밀도, 마하수 항력, 코리올리 힘을 반영한 미분 가속도 산출 (0-Alloc)
        /// 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static (Vector3 dPos, Vector3 dVel) ComputeDerivatives(
            in BallisticState state,
            float mass,
            float referenceArea,
            float dragCoefficient,
            float latitudeRad,
            in Vector3 wind)
        {
            // 1. 공기 상대 속도
            Vector3 vRel = state.Velocity - wind;
            float speedRel = vRel.Length();

            // 2. 1976 미국 표준 대기 간이 밀도 계산 rho(z)
            float altitude = Math.Max(0f, state.Position.Z);
            float rho = 1.225f * MathF.Exp(-altitude / 8435.0f);

            // 3. 공기 항력 가속도 a_drag = -0.5 * rho * A * Cd * vRel * vRel_vec / m
            float dragMagnitude = 0.5f * rho * referenceArea * dragCoefficient * speedRel / mass;
            Vector3 aDrag = -dragMagnitude * vRel;

            // 4. 중력 가속도 (Z축 아래)
            Vector3 aGravity = new Vector3(0, 0, -StandardGravity);

            // 5. 3차원 코리올리 가속도 a_coriolis = -2 * (omega x v)
            // omega = (0, EarthAngularSpeed * cos(lat), EarthAngularSpeed * sin(lat))
            Vector3 omega = new Vector3(
                0f,
                EarthAngularSpeed * MathF.Cos(latitudeRad),
                EarthAngularSpeed * MathF.Sin(latitudeRad)
            );
            Vector3 aCoriolis = -2f * Vector3.Cross(omega, state.Velocity);

            Vector3 totalAcc = aGravity + aDrag + aCoriolis;
            return (state.Velocity, totalAcc);
        }

        /// 

        /// 4차 룬게-쿠타(RK4) 단일 시간 스텝 수치 전진
        /// 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static BallisticState StepRK4(
            in BallisticState current,
            float dt,
            float mass,
            float area,
            float cd,
            float latitudeRad,
            in Vector3 wind)
        {
            var (dp1, dv1) = ComputeDerivatives(current, mass, area, cd, latitudeRad, wind);

            BallisticState s2 = new BallisticState(
                current.Position + dp1 * (dt * 0.5f),
                current.Velocity + dv1 * (dt * 0.5f));
            var (dp2, dv2) = ComputeDerivatives(s2, mass, area, cd, latitudeRad, wind);

            BallisticState s3 = new BallisticState(
                current.Position + dp2 * (dt * 0.5f),
                current.Velocity + dv2 * (dt * 0.5f));
            var (dp3, dv3) = ComputeDerivatives(s3, mass, area, cd, latitudeRad, wind);

            BallisticState s4 = new BallisticState(
                current.Position + dp3 * dt,
                current.Velocity + dv3 * dt);
            var (dp4, dv4) = ComputeDerivatives(s4, mass, area, cd, latitudeRad, wind);

            Vector3 nextPos = current.Position + (dp1 + 2f * dp2 + 2f * dp3 + dp4) * (dt / 6.0f);
            Vector3 nextVel = current.Velocity + (dv1 + 2f * dv2 + 2f * dv3 + dv4) * (dt / 6.0f);

            return new BallisticState(nextPos, nextVel);
        }
    }
}