clc;
clear;
load("20210106150614_01_staring.mat")

% 假设已经加载了数据并有频率数据freq
measured_sigma0 = calculateScatteringCoefficient(amplitude_complex_T1); % 使用T1数据计算实测散射系数

% 假设已经有模型的数据
freq = [1e9, 2e9, 3e9]; % 频率数组示例

%% GIT
GrAng=0.3;
alpha = GrAng;
ThWind=0;
Pol='H';
SeaState=3;
for ifreq = 1 : length(freq)
GTI_SigZHor(ifreq) = GTI_SigmaSea(freq(ifreq),SeaState,Pol,alpha,ThWind);
end
GTI_SigmaHH = GTI_SigZHor;

%% NRL
NRL_SigZHor = NRL_SigmaSea(freq,SeaState,Pol,alpha);
NRL_SigmaHH=NRL_SigZHor;

%% HYB
for ifreq = 1 : length(freq)
HYB_SigZHor(ifreq) = HYB_SigmaSea(freq(ifreq),SeaState,Pol,alpha,ThWind);
end
HYB_SigmaHH=HYB_SigZHor;
%% TSC
for ifreq = 1 : length(freq)
TSC_SigZHor(ifreq) = TSC_SigmaSea(freq(ifreq),SeaState,Pol,alpha,ThWind);
end
TSC_SigmaHH=TSC_SigZHor;
%%
modelSigmas.GTI = GTI_SigmaHH; % GIT模型数据
modelSigmas.NRL = NRL_SigmaHH; % NRL模型数据
modelSigmas.HYB = HYB_SigmaHH; % HYB模型数据
modelSigmas.TSC = TSC_SigmaHH; % TSC模型数据
% 绘制比较图
plotScatteringCoefficients(freq, modelSigmas, measured_sigma0);

function sigma0 = calculateScatteringCoefficient(amplitude_complex)
% 计算每个点的功率
power_spectrum = abs(amplitude_complex).^2; % I^2 + Q^2

% 设置参考功率,例如1毫瓦(1e-3瓦特)
P_ref = 1e-3;

% 计算散射系数(dB)
sigma0 = 10 * log10(power_spectrum / P_ref);  % dB

end

function plotScatteringCoefficients(freq, modelSigmas, measuredSigma)
% freq: 频率数组
% modelSigmas: 结构体,包含不同模型的散射系数数据
% measuredSigma: 实测散射系数数组
figure;
semilogx(freq, modelSigmas.GTI, 'r-', freq, modelSigmas.NRL, 'k-', freq, modelSigmas.HYB, 'b-', freq, modelSigmas.TSC, 'g-', freq, measuredSigma, 'm--');
hold on;
legend('GIT', 'NRL', 'HYB', 'TSC', 'Measured');
xlabel('Frequency (Hz)');
ylabel('Scattering Coefficient (dB)');
title('Comparison of Scattering Coefficient Models and Measurement');
grid on;
end