초장기선 전파간섭계(VLBI) 위상 보정 & 1,000Hz 0-Alloc C# 반시테르트-체르니케 2D 역푸리에 천체 영상 합성 엔진

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

namespace ZeliDesk.Education.StemAstro.VlbiInterferometry
{
    public readonly record struct BaselineVisibility(
        float U_Wavelengths, // 기선 공간주파수 U (파장 단위)
        float V_Wavelengths, // 기선 공간주파수 V (파장 단위)
        Complex RawVisibility, // 측정된 복소 가시도
        float Station1AtmPhaseRad, // 안테나 1 대기 위상 오차
        float Station2AtmPhaseRad  // 안테나 2 대기 위상 오차
    );

    public ref struct SynthesizedPixel
    {
        public float DirectionL; // 천체 각도 좌표 L (rad)
        public float DirectionM; // 천체 각도 좌표 M (rad)
        public float BrightnessJansky; // 합성 밝기 (Jy/beam)
    }

    public static class VlbiSynthesizerCore
    {
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static Complex CalibrateAtmosphericPhase(in BaselineVisibility vis)
        {
            // 1. 대기 요동 위상차 역회전 보정 (0-Alloc 복소수 연산)
            float netPhaseError = vis.Station1AtmPhaseRad - vis.Station2AtmPhaseRad;
            Complex phaseRotator = Complex.FromPolarCoordinates(1.0, -netPhaseError);
            return vis.RawVisibility * phaseRotator;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void SynthesizePixelBatch(
            ReadOnlySpan visibilities,
            float lRad,
            float mRad,
            ref SynthesizedPixel output)
        {
            output.DirectionL = lRad;
            output.DirectionM = mRad;

            float realAccum = 0f;
            float twoPi = 2f * MathF.PI;

            // 2. 반시테르트-체르니케 2D 역푸리에 단일 화소 누적 (Dirty Map)
            for (int i = 0; i < visibilities.Length; i++)
            {
                ref readonly var vis = ref visibilities[i];
                Complex calibrated = CalibrateAtmosphericPhase(in vis);

                float phase = twoPi * (vis.U_Wavelengths * lRad + vis.V_Wavelengths * mRad);
                float cosP = MathF.Cos(phase);
                float sinP = MathF.Sin(phase);

                // e^(i * phase) = cos + i*sin 과의 내적의 실수부
                realAccum += (float)calibrated.Real * cosP - (float)calibrated.Imaginary * sinP;
            }

            output.BrightnessJansky = MathF.Max(0f, realAccum / Math.Max(1, visibilities.Length));
        }
    }
}