function [snr_coverage_2d, azimuth_deg, elevation_deg, coverage_info, fig] = analyze_coverage(RadarParams, target_range_m, target_rcs_dBsm)
    % Coverage 분석: 2D 안테나 패턴 기반 각도별 SNR 계산
    % 
    % 방법:
    % 1) Azimuth 0도, Elevation 0도에서 기준 SNR 계산
    % 2) 다른 각도의 상대 이득 차이로 SNR 계산 (효율적)
    % 3) 입력한 RCS 값 반영
    %
    % 입력:
    %   - RadarParams: 레이더 파라미터 구조체
    %   - target_range_m: 표적 거리 (m), 기본값 100
    %   - target_rcs_dBsm: 표적 RCS (dBsm), 기본값 0
    %
    % 출력:
    %   - snr_coverage_2d: 2D SNR 맵 [elevation x azimuth]
    %   - azimuth_deg: 방위각 벡터
    %   - elevation_deg: 고각 벡터
    %   - coverage_info: 커버리지 정보 구조체
    %   - fig: 생성된 figure 핸들
    
    %% 입력 파라미터 처리
    if nargin < 2
        target_range_m = 100;  % 기본: 100m
    end
    if nargin < 3
        target_rcs_dBsm = 0;   % 기본: 0 dBsm = 1
    end
    
    %% 기본 파라미터 추출
    lambda = RadarParams.Waveform.lambda_c;
    fc = RadarParams.Waveform.fc;
    kb = RadarParams.Basic.kb;
    T0 = RadarParams.Basic.T0;
    
    % 안테나 및 경로 파라미터
    TxPattern = RadarParams.Antenna.TxPattern;
    RxPattern = RadarParams.Antenna.RxPattern;
    NumTx = RadarParams.Antenna.NumTx;
    NumRx = RadarParams.Antenna.NumRx;
    
    % 수신 경로 파라미터
    rxPathGain_dB = RadarParams.Rxpath.rxPathGain_dB;
    system_NF_dB = RadarParams.Rxpath.system_NF_dB;
    
    % RCS (선형값)
    rcs_target = 10^(target_rcs_dBsm / 10);
    
    %% 평균 2D 안테나 패턴 생성 (1D -> 2D 변환 포함)
    % 공통 그리드 설정 (높은 해상도)
    az_common = -90:1:90;
    el_common = -90:1:90;
    [AZ_common_grid, EL_common_grid] = meshgrid(az_common, el_common);
    
    % TX 패턴들을 2D로 변환 및 평균
    avg_tx_gain_2d = zeros(length(el_common), length(az_common));
    for tx = 1:NumTx
        if isfield(TxPattern(tx), 'gain_az_dBi')  % 1D 패턴
            az_angles = TxPattern(tx).az_angles;
            el_angles = TxPattern(tx).el_angles;
            gain_az = TxPattern(tx).gain_az_dBi;
            gain_el = TxPattern(tx).gain_el_dBi;
            
            % 각 패턴의 최대값 저장
            max_gain_az = max(gain_az);
            max_gain_el = max(gain_el);
            max_gain_ref = max(max_gain_az, max_gain_el);
            
            % Normalize (최대값 = 0 dB)
            gain_az_norm = gain_az - max_gain_az;
            gain_el_norm = gain_el - max_gain_el;
            
            % 정규화된 패턴 보간
            gain_az_interp = interp1(az_angles, gain_az_norm, AZ_common_grid, 'linear', 'extrap');
            gain_el_interp = interp1(el_angles, gain_el_norm, EL_common_grid, 'linear', 'extrap');
            
            % 2D로 결합하고 최대 이득 더하기
            tx_gain_2d = gain_az_interp + gain_el_interp + max_gain_ref;
        else  % 2D 패턴
            az_angles = TxPattern(tx).az_angles;
            el_angles = TxPattern(tx).el_angles;
            gain_dBi = TxPattern(tx).gain_dBi;
            [AZ_grid, EL_grid] = ndgrid(az_angles, el_angles);
            F = scatteredInterpolant(AZ_grid(:), EL_grid(:), gain_dBi(:), 'linear', 'nearest');
            tx_gain_2d = F(AZ_common_grid, EL_common_grid);
        end
        avg_tx_gain_2d = avg_tx_gain_2d + tx_gain_2d;
    end
    avg_tx_gain_2d = avg_tx_gain_2d / NumTx;
    
    % RX 패턴들을 2D로 변환 및 평균 (동일한 방식)
    avg_rx_gain_2d = zeros(length(el_common), length(az_common));
    for rx = 1:NumRx
        if isfield(RxPattern(rx), 'gain_az_dBi')  % 1D 패턴
            az_angles = RxPattern(rx).az_angles;
            el_angles = RxPattern(rx).el_angles;
            gain_az = RxPattern(rx).gain_az_dBi;
            gain_el = RxPattern(rx).gain_el_dBi;
            
            % 각 패턴의 최대값 저장
            max_gain_az = max(gain_az);
            max_gain_el = max(gain_el);
            max_gain_ref = max(max_gain_az, max_gain_el);
            
            % Normalize (최대값 = 0 dB)
            gain_az_norm = gain_az - max_gain_az;
            gain_el_norm = gain_el - max_gain_el;
            
            % 정규화된 패턴 보간
            gain_az_interp = interp1(az_angles, gain_az_norm, AZ_common_grid, 'linear', 'extrap');
            gain_el_interp = interp1(el_angles, gain_el_norm, EL_common_grid, 'linear', 'extrap');
            
            % 2D로 결합하고 최대 이득 더하기
            rx_gain_2d = gain_az_interp + gain_el_interp + max_gain_ref;
        else  % 2D 패턴
            az_angles = RxPattern(rx).az_angles;
            el_angles = RxPattern(rx).el_angles;
            gain_dBi = RxPattern(rx).gain_dBi;
            [AZ_grid, EL_grid] = ndgrid(az_angles, el_angles);
            F = scatteredInterpolant(AZ_grid(:), EL_grid(:), gain_dBi(:), 'linear', 'nearest');
            rx_gain_2d = F(AZ_common_grid, EL_common_grid);
        end
        avg_rx_gain_2d = avg_rx_gain_2d + rx_gain_2d;
    end
    avg_rx_gain_2d = avg_rx_gain_2d / NumRx;
    
    %% TX 전력 및 공통 파라미터 계산
    % TX 전력
    pa_profile_freq = RadarParams.RFOutput.PA_Profile.freqs;
    pa_profile_power_dbm = RadarParams.RFOutput.PA_Profile.power_dBm;
    ptx_dbm = interp1(pa_profile_freq, pa_profile_power_dbm, fc, 'linear', 'extrap');
    ptx_w = 10^((ptx_dbm - 30) / 10);
    
    % 경로 손실 및 잡음 계산 (모든 각도에서 공통)
    path_loss_factor = (4 * pi * target_range_m)^2;
    rxGain_linear = 10^(rxPathGain_dB / 10);
    noise_power_w = kb * T0 * RadarParams.Waveform.fs_adc;
    system_NF_linear = 10^(system_NF_dB / 10);
    total_noise_power_w = noise_power_w * system_NF_linear * rxGain_linear;  % RX 이득 포함
    
    %% 각도별 SNR 계산 (2D 안테나 패턴 직접 사용 - 벡터화 연산)
    % X축: Azimuth, Y축: Elevation, Z축: SNR
    azimuth_deg = -90:1:90;   % X축: -90 ~ 90도, 1도 단위
    elevation_deg = -90:1:90;  % Y축: -90 ~ 90도, 1도 단위
    
    % az_common과 elevation_deg가 동일한 그리드이므로 인덱스 매칭
    % 안테나 이득 추출 (dB)
    g_tx_db_grid = avg_tx_gain_2d;  % [num_el x num_az]
    g_rx_db_grid = avg_rx_gain_2d;  % [num_el x num_az]
    
    % 이득을 선형으로 변환
    g_tx_linear_grid = 10.^(g_tx_db_grid / 10);
    g_rx_linear_grid = 10.^(g_rx_db_grid / 10);

    % 신호처리 이득
    SP_gain_rngFFT = RadarParams.Waveform.Timing.AdcSampTime * RadarParams.Waveform.fs_adc;
    SP_gain_dopFFT = RadarParams.Waveform.NumChirps;
    SP_loss_rngwin = 10^(RadarParams.SP.RDM.window_metrics_range.snr_loss_dB / 10);
    SP_loss_dopwin = 10^(RadarParams.SP.RDM.window_metrics_doppler.snr_loss_dB / 10);
    SP_loss_rng_straddle = 10^(RadarParams.SP.RDM.window_metrics_range.scalloping_loss_dB / 10);
    SP_loss_dop_straddle = 10^(RadarParams.SP.RDM.window_metrics_doppler.scalloping_loss_dB / 10);
    SP_total = SP_gain_rngFFT * SP_gain_dopFFT / (SP_loss_rngwin * SP_loss_dopwin * SP_loss_rng_straddle * SP_loss_dop_straddle);

    % Secondary Surface Loss
    secondary_loss = 10^(RadarParams.Antenna.SecondarySurfaceLoss_dB / 10);

    fprintf('SP gain rng FFT: %.2f, SP gain doppler FFT: %.2f\n', SP_gain_rngFFT, SP_gain_dopFFT);
    fprintf('SP loss range window: %.2f dB, SP loss doppler window: %.2f dB\n', RadarParams.SP.RDM.window_metrics_range.snr_loss_dB, RadarParams.SP.RDM.window_metrics_doppler.snr_loss_dB);
    fprintf('SP loss range straddle: %.2f dB, SP loss doppler straddle: %.2f dB\n', RadarParams.SP.RDM.window_metrics_range.scalloping_loss_dB, RadarParams.SP.RDM.window_metrics_doppler.scalloping_loss_dB);
    fprintf('SP loss doppler window: %.2f dB\n', RadarParams.SP.RDM.window_metrics_doppler.snr_loss_dB);
    fprintf('SP loss doppler straddle: %.2f dB\n', RadarParams.SP.RDM.window_metrics_doppler.scalloping_loss_dB);
    fprintf('Secondary surface loss: %.2f dB\n', RadarParams.Antenna.SecondarySurfaceLoss_dB);
    
    % SNR 계산: 벡터화 연산 (스칼라 항 × 2D 배열)
    p_rx_w = (ptx_w * (lambda^2) * rcs_target) / ((path_loss_factor^2) * (4 * pi)) * SP_total / secondary_loss;
    p_rx_w = p_rx_w * g_tx_linear_grid .* g_rx_linear_grid; % 연산량 최적화를 위해 안테나 이득을 나중에 곱함
    p_rx_after_rxgain_w = p_rx_w * rxGain_linear;
    snr_linear = p_rx_after_rxgain_w / total_noise_power_w;
    snr_coverage_2d = 10 * log10(max(snr_linear, eps));
    
    %% 1D CUT 추출 (elevation=0에서의 Azimuth cut, azimuth=0에서의 Elevation cut)
    [~, el_idx_zero] = min(abs(elevation_deg - 0));
    snr_az_cut = snr_coverage_2d(el_idx_zero, :);  % elevation=0
    
    [~, az_idx_zero] = min(abs(azimuth_deg - 0));
    snr_el_cut = snr_coverage_2d(:, az_idx_zero);  % azimuth=0
    
    %% Cartesian 좌표계 변환 (Range, Azimuth, Elevation -> x, y, z)
    [AZ_deg_grid, EL_deg_grid] = meshgrid(azimuth_deg, elevation_deg);
    x_grid_m = target_range_m .* cosd(EL_deg_grid) .* cosd(AZ_deg_grid);
    y_grid_m = target_range_m .* cosd(EL_deg_grid) .* sind(AZ_deg_grid);
    z_grid_m = target_range_m .* sind(EL_deg_grid);

    %% Coverage 정보 통계
    snr_2d_vec = snr_coverage_2d(:);
    
    coverage_info.azimuth_deg = azimuth_deg;
    coverage_info.elevation_deg = elevation_deg;
    coverage_info.snr_coverage_2d = snr_coverage_2d;
    coverage_info.snr_az_cut = snr_az_cut;
    coverage_info.snr_el_cut = snr_el_cut;
    coverage_info.target_range_m = target_range_m;
    coverage_info.target_rcs_dBsm = target_rcs_dBsm;
    coverage_info.x_grid_m = x_grid_m;
    coverage_info.y_grid_m = y_grid_m;
    coverage_info.z_grid_m = z_grid_m;
    coverage_info.mean_snr = mean(snr_2d_vec);
    coverage_info.max_snr = max(snr_2d_vec);
    coverage_info.min_snr = min(snr_2d_vec);
    [~, max_idx_2d] = max(snr_2d_vec);
    [max_el_idx_2d, max_az_idx_2d] = ind2sub(size(snr_coverage_2d), max_idx_2d);
    coverage_info.max_snr_azimuth = azimuth_deg(max_az_idx_2d);
    coverage_info.max_snr_elevation = elevation_deg(max_el_idx_2d);
    coverage_info.max_snr_x_m = x_grid_m(max_el_idx_2d, max_az_idx_2d);
    coverage_info.max_snr_y_m = y_grid_m(max_el_idx_2d, max_az_idx_2d);
    coverage_info.max_snr_z_m = z_grid_m(max_el_idx_2d, max_az_idx_2d);
    
    %% 3D Coverage 시각화
    fig = figure('Name', 'Coverage Analysis - 3D');
    set(fig, 'Position', [100, 100, 1600, 900]);
    
    % 1D cut SNR 추출
    [~, el_idx_zero] = min(abs(elevation_deg - 0));
    [~, az_idx_zero] = min(abs(azimuth_deg - 0));
    snr_az_cut = snr_coverage_2d(el_idx_zero, :);  % elevation=0 cut
    snr_el_cut = snr_coverage_2d(:, az_idx_zero);  % azimuth=0 cut

    % Subplot 1: x-y top view (color = SNR at elevation=0 cut)
    ax_topview = subplot(2, 2, 3);
    x_top_m = target_range_m .* cosd(azimuth_deg);
    y_top_m = target_range_m .* sind(azimuth_deg);
    scatter(x_top_m, y_top_m, 45, snr_az_cut, 'filled');
    axis equal;
    grid on;
    cb = colorbar;
    ylabel(cb, 'SNR (dB)');
    xlabel('x (m)');
    ylabel('y (m)');
    title('Top View (Azimuth)');
    
    % Subplot 2: 3D Surface (geometry = x,y,z, color = SNR) - spans full width at top
    ax_3d = subplot(2, 2, 1:2);
    surf(x_grid_m, y_grid_m, z_grid_m, snr_coverage_2d, 'EdgeColor', 'none');
    colormap(jet);
    cb_3d = colorbar;
    ylabel(cb_3d, 'SNR (dB)');
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('z (m)');
    title('3D Coverage in Cartesian Coordinates');
    view(45, 30);
    grid on;
    hold on;
    plot3(coverage_info.max_snr_x_m, coverage_info.max_snr_y_m, coverage_info.max_snr_z_m, 'r*', 'MarkerSize', 20, 'LineWidth', 2);
    hold off;
    
    % Subplot 3: x-z side view (azimuth=0, color=SNR at azimuth=0 cut)
    ax_sideview = subplot(2, 2, 4);
    x_side_m = target_range_m .* cosd(elevation_deg);
    z_side_m = target_range_m .* sind(elevation_deg);
    scatter(x_side_m, z_side_m, 45, snr_el_cut, 'filled');
    grid on;
    cb_side = colorbar;
    ylabel(cb_side, 'SNR (dB)');
    xlabel('x (m)');
    ylabel('z (m)');
    title('Side View (Elevation)');

    % 전체 타이틀
    sgtitle(sprintf('3D Coverage (x,y,z) | Max SNR: %.2f dB at (%.2f, %.2f, %.2f) m', ...
                    coverage_info.max_snr, coverage_info.max_snr_x_m, coverage_info.max_snr_y_m, coverage_info.max_snr_z_m), 'FontSize', 12);
    
    % 단일 Data Cursor Mode 설정 (모든 축에 적용)
    dcm_obj = datacursormode(fig);
    dcm_obj.Enable = 'on';
    set(dcm_obj, 'UpdateFcn', {@unified_datatip_callback, ...
        ax_topview, ax_3d, ax_sideview, ...
        x_top_m, y_top_m, snr_az_cut, azimuth_deg, ...
        x_grid_m, y_grid_m, z_grid_m, snr_coverage_2d, elevation_deg, ...
        x_side_m, z_side_m, snr_el_cut});
    dcm_obj.SnapToDataVertex = 'on';
    
    % Help text
    fprintf('\n>>> Data Cursor Mode ENABLED <<<\n');
    fprintf('Instructions:\n');
    fprintf('  1. Left-click on any plot to inspect points\n');
    fprintf('  2. Top View shows: x, y, Azimuth, SNR\n');
    fprintf('  3. 3D Surface shows: x, y, z, Azimuth, Elevation, SNR\n');
    fprintf('  4. Side View shows: x, z, Elevation, SNR\n');
    fprintf('  5. Press Escape or click "Disable Data Cursor" in Figure Tools to deactivate\n\n');
    
    %% 통계 출력
    fprintf('\n===== 3D COVERAGE ANALYSIS REPORT (Cartesian x,y,z) =====\n');
    fprintf('Target Configuration:\n');
    fprintf('  - Range: %.1f m\n', target_range_m);
    fprintf('  - RCS: %.1f dBsm\n', target_rcs_dBsm);
    fprintf('  - Analysis Angles: Azimuth ±90 deg, Elevation ±90 deg\n');
    fprintf('\nSNR Statistics (Full Grid):\n');
    fprintf('  - Mean SNR: %.2f dB\n', coverage_info.mean_snr);
    fprintf('  - Max SNR: %.2f dB @ (x,y,z)=(%.2f, %.2f, %.2f) m\n', coverage_info.max_snr, coverage_info.max_snr_x_m, coverage_info.max_snr_y_m, coverage_info.max_snr_z_m);
    fprintf('  - Min SNR: %.2f dB\n', coverage_info.min_snr);
    fprintf('  - SNR Range: %.2f dB\n', coverage_info.max_snr - coverage_info.min_snr);
    fprintf('System Parameters:\n');
    fprintf('  - TX Power: %.2f dBm\n', ptx_dbm);
    fprintf('  - RX Gain: %.1f dB\n', rxPathGain_dB);
    fprintf('  - System NF: %.1f dB\n', system_NF_dB);
    fprintf('=============================================\n\n');

end

function output_txt = unified_datatip_callback(~, event, ...
    ax_topview, ax_3d, ax_sideview, ...
    x_top_m, y_top_m, snr_az_cut, azimuth_deg, ...
    x_grid_m, y_grid_m, z_grid_m, snr_coverage_2d, elevation_deg, ...
    x_side_m, z_side_m, snr_el_cut)
    % Unified datatip callback for all subplots
    
    try
        pos = event.Position;
        current_ax = event.Target.Parent;  % Get current axis
        
        % Determine which subplot was clicked
        if isequal(current_ax, ax_3d)
            % 3D Plot
            x_val = pos(1);
            y_val = pos(2);
            z_val = pos(3);
            
            dist2 = (x_grid_m - x_val).^2 + (y_grid_m - y_val).^2 + (z_grid_m - z_val).^2;
            [~, nearest_idx] = min(dist2(:));
            [nearest_el_idx, nearest_az_idx] = ind2sub(size(dist2), nearest_idx);
            snr_val = snr_coverage_2d(nearest_el_idx, nearest_az_idx);
            az_val = azimuth_deg(nearest_az_idx);
            el_val = elevation_deg(nearest_el_idx);
            
            output_txt = {
                ['x: ' num2str(x_val, '%.3f') ' m']
                ['y: ' num2str(y_val, '%.3f') ' m']
                ['z: ' num2str(z_val, '%.3f') ' m']
                ['Azimuth: ' num2str(az_val, '%.2f') ' deg']
                ['Elevation: ' num2str(el_val, '%.2f') ' deg']
                ['SNR[dB]: ' num2str(snr_val, '%.3f')]
            };
        elseif isequal(current_ax, ax_topview)
            % Top View - x, y, Azimuth, SNR
            x_val = pos(1);
            y_val = pos(2);
            
            [~, nearest_az_idx] = min(abs(x_top_m - x_val).^2 + abs(y_top_m - y_val).^2);
            az_val = azimuth_deg(nearest_az_idx);
            snr_val = snr_az_cut(nearest_az_idx);
            
            output_txt = {
                ['x: ' num2str(x_val, '%.3f') ' m']
                ['y: ' num2str(y_val, '%.3f') ' m']
                ['Azimuth: ' num2str(az_val, '%.2f') ' deg']
                ['SNR[dB]: ' num2str(snr_val, '%.3f')]
            };
        elseif isequal(current_ax, ax_sideview)
            % Side View - x, z, Elevation, SNR
            x_val = pos(1);
            z_val = pos(2);
            
            [~, nearest_el_idx] = min(abs(x_side_m - x_val).^2 + abs(z_side_m - z_val).^2);
            el_val = elevation_deg(nearest_el_idx);
            snr_val = snr_el_cut(nearest_el_idx);
            
            output_txt = {
                ['x: ' num2str(x_val, '%.3f') ' m']
                ['z: ' num2str(z_val, '%.3f') ' m']
                ['Elevation: ' num2str(el_val, '%.2f') ' deg']
                ['SNR[dB]: ' num2str(snr_val, '%.3f')]
            };
        else
            output_txt = 'Unknown plot';
        end
    catch ME
        output_txt = ['Error: ' ME.message];
    end
end