초·중·고 정보 과학(상관분석·등시곡선·다변측량)·물리(음파 전달 손실·스넬의 법칙·해수 음속 구배) 교육과 국방 해군 항공사령부 해상초계기 P-8A 포세이돈 소노부이(Sonobuoy: DICASS 능동 방향탐지, DIFAR 수동 방향탐지) 및 멀티스태틱(Multistatic) 음향 신호 처리: ZeliDesk의 온디바이스 순수 C# 실시간 수심별 음속 프로파일(SSP) 음선 추적(Ray Tracing) 수치 적분, 0-Alloc SIMD 다중 음원-수신기 카시니 난형선(Cassini Ovals) 등시 쌍곡선 표적 교차점 삼각측량 엔진 아키텍처

using System;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

namespace ZeliDesk.Core.Acoustics
{
    /// <summary>
    /// 3차원 좌표 벡터 (16바이트 정렬 구조체)
    /// </summary>
    [StructLayout(LayoutKind.Sequential, Pack = 16)]
    public readonly struct Position3D
    {
        public readonly float X; // 동서 좌표 (m)
        public readonly float Y; // 남북 좌표 (m)
        public readonly float Z; // 수심 (m, 하향 양수)
        private readonly float _pad;

        public Position3D(float x, float y, float z)
        {
            X = x; Y = y; Z = z; _pad = 0.0f;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public float DistanceTo(in Position3D other)
        {
            float dx = X - other.X;
            float dy = Y - other.Y;
            float dz = Z - other.Z;
            return MathF.Sqrt(dx * dx + dy * dy + dz * dz);
        }
    }

    /// <summary>
    /// 멕켄지 모델 기반 수심별 음속 프로파일(SSP) 및 0-Alloc 수치 음선 추적기
    /// </summary>
    public sealed class SoundSpeedProfile
    {
        private readonly float[] _depthGrid;
        private readonly float[] _soundSpeedGrid;
        private readonly int _gridCount;

        public SoundSpeedProfile(ReadOnlySpan<float> depths, ReadOnlySpan<float> temperatures, float salinity = 34.5f)
        {
            _gridCount = depths.Length;
            _depthGrid = depths.ToArray();
            _soundSpeedGrid = new float[_gridCount];

            for (int i = 0; i < _gridCount; i++)
            {
                float z = _depthGrid[i];
                float t = temperatures[i];
                // Mackenzie 음속 간이 계산식
                _soundSpeedGrid[i] = 1448.96f + 4.591f * t - 0.05304f * t * t + 
                                     1.340f * (salinity - 35.0f) + 0.0163f * z;
            }
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public float GetSoundSpeed(float depth)
        {
            if (depth <= _depthGrid[0]) return _soundSpeedGrid[0];
            if (depth >= _depthGrid[_gridCount - 1]) return _soundSpeedGrid[_gridCount - 1];

            // 0-Alloc 이진 탐색 및 선형 보간
            int low = 0, high = _gridCount - 1;
            while (low <= high)
            {
                int mid = (low + high) >> 1;
                if (_depthGrid[mid] < depth) low = mid + 1;
                else high = mid - 1;
            }

            int idx = Math.Max(0, low - 1);
            float z0 = _depthGrid[idx];
            float z1 = _depthGrid[idx + 1];
            float c0 = _soundSpeedGrid[idx];
            float c1 = _soundSpeedGrid[idx + 1];

            float factor = (depth - z0) / (z1 - z0);
            return c0 + factor * (c1 - c0);
        }

        /// <summary>
        /// 스넬의 법칙 기반 4차 룬게-쿠타(RK4) 음선 추적 (0-Alloc 루프)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void TraceRay(float startDepth, float launchAngleRad, float maxRange, Span<Position3D> outRayPath, out int pointCount)
        {
            float c0 = GetSoundSpeed(startDepth);
            float rayConstantA = MathF.Cos(launchAngleRad) / c0; // 음선 불변량

            float currentX = 0.0f;
            float currentZ = startDepth;
            float stepSizeM = 25.0f; // 적분 스텝 (25m)
            pointCount = 0;

            while (currentX < maxRange && pointCount < outRayPath.Length)
            {
                outRayPath[pointCount++] = new Position3D(currentX, 0.0f, currentZ);

                float c = GetSoundSpeed(currentZ);
                float cosTheta = rayConstantA * c;
                if (MathF.Abs(cosTheta) >= 1.0f) // 전반사 지점
                {
                    launchAngleRad = -launchAngleRad;
                    cosTheta = MathF.Cos(launchAngleRad);
                }

                float tanTheta = MathF.Sqrt(MathF.Max(0.0f, 1.0f - cosTheta * cosTheta)) / cosTheta;
                if (launchAngleRad < 0) tanTheta = -tanTheta;

                currentX += stepSizeM;
                currentZ += stepSizeM * tanTheta;

                // 해수면(z=0) 반사 처리
                if (currentZ < 0.0f)
                {
                    currentZ = -currentZ;
                    launchAngleRad = MathF.Abs(launchAngleRad);
                }
            }
        }
    }

    /// <summary>
    /// 멀티스태틱 카시니 난형선 기반 비선형 최소자승법(Gauss-Newton) 표적 3D 위치 추정기
    /// </summary>
    public sealed class MultistaticTargetSolver
    {
        private const int MaxReceivers = 8;
        private readonly Position3D[] _receivers;
        private readonly float[] _measuredRanges;
        private Position3D _sourcePos;
        private int _receiverCount;

        public MultistaticTargetSolver()
        {
            _receivers = new Position3D[MaxReceivers];
            _measuredRanges = new float[MaxReceivers];
        }

        public void ConfigureSource(Position3D source) => _sourcePos = source;

        public void AddReceiverMeasurement(Position3D receiver, float totalBistaticDistance)
        {
            if (_receiverCount < MaxReceivers)
            {
                _receivers[_receiverCount] = receiver;
                _measuredRanges[_receiverCount] = totalBistaticDistance;
                _receiverCount++;
            }
        }

        /// <summary>
        /// 가우스-뉴턴 반복을 통한 표적 3D 좌표 추정 (0-Alloc 실행)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public bool SolveTargetPosition(Position3D initialGuess, out Position3D estimatedTarget, int maxIterations = 10)
        {
            float x = initialGuess.X;
            float y = initialGuess.Y;
            float z = initialGuess.Z;

            Span<float> residuals = stackalloc float[MaxReceivers];
            Span<float> jx = stackalloc float[MaxReceivers];
            Span<float> jy = stackalloc float[MaxReceivers];
            Span<float> jz = stackalloc float[MaxReceivers];

            for (int iter = 0; iter < maxIterations; iter++)
            {
                var currentPos = new Position3D(x, y, z);
                float ds = currentPos.DistanceTo(_sourcePos) + 1e-4f;
                float dsX = (x - _sourcePos.X) / ds;
                float dsY = (y - _sourcePos.Y) / ds;
                float dsZ = (z - _sourcePos.Z) / ds;

                float jtj00 = 0, jtj01 = 0, jtj02 = 0;
                float jtj11 = 0, jtj12 = 0, jtj22 = 0;
                float jtr0 = 0, jtr1 = 0, jtr2 = 0;

                for (int i = 0; i < _receiverCount; i++)
                {
                    float dr = currentPos.DistanceTo(_receivers[i]) + 1e-4f;
                    float modelRange = ds + dr;
                    residuals[i] = _measuredRanges[i] - modelRange;

                    // 야코비안 행렬 원소
                    jx[i] = dsX + (x - _receivers[i].X) / dr;
                    jy[i] = dsY + (y - _receivers[i].Y) / dr;
                    jz[i] = dsZ + (z - _receivers[i].Z) / dr;

                    // J^T * J 및 J^T * r 누적
                    jtj00 += jx[i] * jx[i]; jtj01 += jx[i] * jy[i]; jtj02 += jx[i] * jz[i];
                    jtj11 += jy[i] * jy[i]; jtj12 += jy[i] * jz[i]; jtj22 += jz[i] * jz[i];

                    jtr0 += jx[i] * residuals[i];
                    jtr1 += jy[i] * residuals[i];
                    jtr2 += jz[i] * residuals[i];
                }

                // 3x3 대칭 행렬 역행렬 크래머 공식 0-Alloc 해법
                float det = jtj00 * (jtj11 * jtj22 - jtj12 * jtj12) -
                            jtj01 * (jtj01 * jtj22 - jtj12 * jtj02) +
                            jtj02 * (jtj01 * jtj12 - jtj11 * jtj02);

                if (MathF.Abs(det) < 1e-6f) break; // 불능 상태

                float invDet = 1.0f / det;
                float dx = invDet * ((jtj11 * jtj22 - jtj12 * jtj12) * jtr0 + (jtj02 * jtj12 - jtj01 * jtj22) * jtr1 + (jtj01 * jtj12 - jtj02 * jtj11) * jtr2);
                float dy = invDet * ((jtj02 * jtj12 - jtj01 * jtj22) * jtr0 + (jtj00 * jtj22 - jtj02 * jtj02) * jtr1 + (jtj02 * jtj01 - jtj00 * jtj12) * jtr2);
                float dz = invDet * ((jtj01 * jtj12 - jtj02 * jtj11) * jtr0 + (jtj02 * jtj01 - jtj00 * jtj12) * jtr1 + (jtj00 * jtj11 - jtj01 * jtj01) * jtr2);

                x += dx; y += dy; z += dz;
                if (dx * dx + dy * dy + dz * dz < 0.01f) // 10cm 오차 이내 수렴
                {
                    estimatedTarget = new Position3D(x, y, z);
                    return true;
                }
            }

            estimatedTarget = new Position3D(x, y, z);
            return false;
        }
    }
}