지구 액체 외핵 지오다이나모(Geodynamo) 쌍극자 역전 & 1,000Hz 0-Alloc C# 코리올리-로런츠(MAC) 힘 평형 및 구면조화함수 지자기 영년변화 항법 서보 엔진

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

namespace ZeliDesk.Stem.Geomagnetism
{
    /// 

    /// 1,000Hz 지구 외핵 지오다이나모 지자기 구면조화 계수 추출 및 SAA 항법 보정 서보
    /// GC 무할당 스택 메모리 0-Alloc 구현체
    /// 
    public ref struct GeodynamoSecularVariationServo
    {
        // 지자기 물리 상수
        private const float EarthMeanRadiusKm = 6371.0f;
        private const float CoreMantleBoundaryRadiusKm = 3480.0f;
        private const float NominalDipoleMomentNt = 29500.0f; // g10 공칭치 (-29,500 nT)
        private const float SaaThresholdTotalFieldNt = 22000.0f; // SAA 이상대 임계치

        // 서보 제어 상태
        public float DipoleGaussG10Nt;
        public float DipoleGaussG11Nt;
        public float DipoleGaussH11Nt;
        public float EstimatedDipoleTiltDeg;
        public float LocalFieldTotalNt;
        public float NavHeadingCorrectionDeg;
        public bool IsInSaaAnomalyZone;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void Initialize()
        {
            DipoleGaussG10Nt = -29450.0f;
            DipoleGaussG11Nt = -1500.0f;
            DipoleGaussH11Nt = 4800.0f;
            EstimatedDipoleTiltDeg = 11.2f;
            LocalFieldTotalNt = 45000.0f;
            NavHeadingCorrectionDeg = 0.0f;
            IsInSaaAnomalyZone = false;
        }

        /// 

        /// 1ms 제어 주기: 위성/항공기 3축 자력계 벡터 및 위경도로부터 가우스 쌍극자 모멘트 및 헤딩 보정각 산출
        /// 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void StepGeomagneticControl(
            float rawBxNt, 
            float rawByNt, 
            float rawBzNt, 
            float latitudeDeg, 
            float longitudeDeg,
            out float currentTotalFieldNt, 
            out float compassDeclinationDeg, 
            out float navHeadingOffsetDeg)
        {
            // 1. 총 자속 밀도 크기 계산 (Zero-Alloc 스택 연산)
            LocalFieldTotalNt = MathF.Sqrt(rawBxNt * rawBxNt + rawByNt * rawByNt + rawBzNt * rawBzNt);
            currentTotalFieldNt = LocalFieldTotalNt;

            // 2. 남대서양 변칙대(SAA) 진입 판별 (총 자계 < 22,000 nT)
            IsInSaaAnomalyZone = LocalFieldTotalNt < SaaThresholdTotalFieldNt;

            // 3. 1차 구면조화함수 가우스 쌍극자 계수 역산 (위경도 기하 투영)
            float latRad = latitudeDeg * (MathF.PI / 180.0f);
            float lonRad = longitudeDeg * (MathF.PI / 180.0f);
            float sinLat = MathF.Sin(latRad);
            float cosLat = MathF.Cos(latRad);

            // Bz = 2 * g10 * sin(lat) 근사
            if (MathF.Abs(sinLat) > 0.1f)
            {
                float estG10 = (rawBzNt * 0.5f) / sinLat;
                DipoleGaussG10Nt = DipoleGaussG10Nt * 0.999f + estG10 * 0.001f;
            }

            // 쌍극자 경사각(Dipole Tilt) 계산: tan(theta) = sqrt(g11^2 + h11^2) / |g10|
            float equatorialDipole = MathF.Sqrt(DipoleGaussG11Nt * DipoleGaussG11Nt + DipoleGaussH11Nt * DipoleGaussH11Nt);
            EstimatedDipoleTiltDeg = MathF.Atan2(equatorialDipole, MathF.Abs(DipoleGaussG10Nt)) * (180.0f / MathF.PI);

            // 4. 자편각(Magnetic Declination D) 계산: D = atan2(By, Bx)
            compassDeclinationDeg = MathF.Atan2(rawByNt, rawBxNt) * (180.0f / MathF.PI);

            // 5. 1,000Hz 관성 항법 자이로 보정 오프셋 생성
            if (IsInSaaAnomalyZone)
            {
                // SAA 지역에서는 자력계 가중치를 줄이고 자이로 드리프트 모델 보정 강화
                NavHeadingCorrectionDeg = compassDeclinationDeg * 0.20f;
            }
            else
            {
                // 일반 지역: 자편각 완전 보상
                NavHeadingCorrectionDeg = compassDeclinationDeg;
            }

            navHeadingOffsetDeg = NavHeadingCorrectionDeg;
        }
    }
}