% 仿真入口，确定 N 和 s 后，通过遍历 eta_1 和 eta_2 来计算出 P

% 设置用于遍历 eta 的参数
eta_1_step = 1e-2;
eta_2_step = 1e-8;

eta_1_start = eta_1_step;
eta_2_start = eta_2_step;

eta_1_end = (sqrt(5) - 1) / 2; % 大于这个值时，eta_1 ^2 + eta_1 必定会大于 1
eta_2_end = 1e-3;

eta_1 = eta_1_start:eta_1_step:eta_1_end;
eta_2 = eta_2_start:eta_2_step:eta_2_end;

eta_1_N = length(eta_1);
eta_2_N = length(eta_2);

N = 1e14;
s = 10;
epsilon = 1e-5;
C_1 = 5576;
C_2 = 10.66;

result = zeros(eta_1_N, eta_2_N);
metric_1 = zeros(eta_1_N, 1);
metric_2 = zeros(eta_2_N, 1);

for i = 1:eta_1_N
    metric_1(i) = (C_1 * s * (log(4*s))^2 * log(8*N)) / (eta_1(i) ^ 2);
end

for j = 1:eta_2_N
    metric_2(j) = (C_2 * s * log(epsilon^-1)) / (eta_2(j) ^ 2);
end


sigmas = zeros(eta_1_N, 1 * eta_2_N);
results = zeros(eta_1_N, 1 * eta_2_N);
k = 0;
d = 1e2;
for i = 1:eta_1_N
    for j = 1:eta_2_N
        sigma = eta_1(i) * (1 + eta_1(i)) + eta_2(j);
        if sigma > 1
            continue
        end
        sigmas(i, j) = sigma;
        k = k + 1;
        m1 = metric_1(i);
        m2 = metric_2(j);

        m = solve_test(m1);
        if m1 < m2
            m = max(m, m2);
        end

        if m > N
            results(i, j) = 0;
        else
            results(i, j) = N / m;
        end
    end
end

figure;
% semilogy(sigmas, results);
h = heatmap(results);
h.GridVisible = false;
ax = gca;
xn = length(ax.XDisplayLabels);
yn = length(ax.YDisplayLabels);

for i = 1:length(ax.XDisplayLabels)
    % if rem(i, rem(xn, 10)) ~= 0
    %     ax.XDisplayLabels(i) = {nan};
    % end
end

for i = 1:length(ax.YDisplayLabels)
    % if rem(i, rem(yn, 10)) ~= 0
    %     ax.YDisplayLabels(i) = {nan};
    % end
end


% heatmap(result);