%% Jacobi矩阵和稳定性分析
function stability_analysis(params_for_condition,condition_label)
    
    fprintf('\n开始稳定性分析...\n');
    
    % 设置参数（使用条件4的参数）
    % params = struct();
    % params.k1 = 1;
    % params.k2 = 1;
    % params.C2 = 5;
    % params.C3 = 100;
    % params.D = 80;
    % params.p1 = 200;
    % params.C = 40;
    % params.m = 50;
    % params.a = 0.5;
    % params.n = 50;
    params=params_for_condition;
    
    % 定义所有平衡点
    equilibrium_points = [
        0, 0, 0;  % E1
        0, 0, 1;  % E2
        0, 1, 0;  % E3
        0, 1, 1;  % E4
        1, 0, 0;  % E5
        1, 0, 1;  % E6
        1, 1, 0;  % E7
        1, 1, 1   % E8
    ];
    
    point_names = {'E1(0,0,0)', 'E2(0,0,1)', 'E3(0,1,0)', 'E4(0,1,1)', ...
                   'E5(1,0,0)', 'E6(1,0,1)', 'E7(1,1,0)', 'E8(1,1,1)'};
    
    % 分析每个平衡点的稳定性
    fprintf('\n平衡点稳定性分析：\n');
    fprintf('%-15s %-20s %-20s %-20s %-15s\n', '平衡点', 'λ1', 'λ2', 'λ3', '稳定性');
    fprintf('%s\n', repmat('-', 1, 90));
    
    stable_points = [];
    
    for i = 1:size(equilibrium_points, 1)
        point = equilibrium_points(i, :);
        
        % 计算Jacobi矩阵
        J = compute_jacobian(point, params);
        
        % 计算特征值
        eigenvalues = eig(J);
        
        % 判断稳定性
        if all(real(eigenvalues) < 0)
            stability = 'ESS (稳定)';
            stable_points = [stable_points; i];
        elseif all(real(eigenvalues) <= 0)
            stability = '中性稳定';
        else
            stability = '不稳定';
        end
        
        % 输出结果
        fprintf('%-15s λ1=%-15.4f λ2=%-15.4f λ3=%-15.4f %-15s\n', ...
                point_names{i}, ...
                real(eigenvalues(1)), ...
                real(eigenvalues(2)), ...
                real(eigenvalues(3)), ...
                stability);
    end
    
    % 可视化稳定性分析
    visualize_stability(equilibrium_points, point_names, params, stable_points);
    
    % 检验条件
    check_stability_conditions(params);
end

function J = compute_jacobian(point, params)
    % 计算给定平衡点的Jacobi矩阵
    x = point(1);
    y = point(2);
    z = point(3);
    
    % 提取参数
    k1 = params.k1;
    k2 = params.k2;
    C2 = params.C2;
    C3 = params.C3;
    D = params.D;
    p1 = params.p1;
    C = params.C;
    m = params.m;
    a = params.a;
    n = params.n;
    
    % Jacobi矩阵元素
    J = zeros(3, 3);
    
    % ?F(x)/?x
    J(1, 1) = (1 - 2*x) * ((C2 * D * y * z) + D - D * y * z - C3 + C2 * z * y);
    
    % ?F(x)/?y
    J(1, 2) = x * (1 - x) * (C2 * D * z - D * z + C2 * z);
    
    % ?F(x)/?z
    J(1, 3) = x * (1 - x) * (C2 * D * y - D * y + C2 * y);
    
    % ?F(y)/?x
    J(2, 1) = y * (1 - y) * m;
    
    % ?F(y)/?y
    J(2, 2) = (1 - 2*y) * (k1 * p1 * z + C * z + m * x + a * n - n - C);
    
    % ?F(y)/?z
    J(2, 3) = y * (1 - y) * (k1 * p1 + C);
    
    % ?F(z)/?x
    J(3, 1) = z * (1 - z) * m;
    
    % ?F(z)/?y
    J(3, 2) = z * (1 - z) * (C + k2 * p1);
    
    % ?F(z)/?z
    J(3, 3) = (1 - 2*z) * (a * n + C * y + m * x + k2 * p1 * y - n - C);
end

function visualize_stability(equilibrium_points, point_names, params, stable_points)
    % 可视化稳定性分析结果
    figure('Position', [100, 100, 1200, 800]);
    
    % 3D相空间中的平衡点
    subplot(2, 2, [1, 3]);
    hold on;
    
    % 绘制所有平衡点
    for i = 1:size(equilibrium_points, 1)
        point = equilibrium_points(i, :);
        if ismember(i, stable_points)
            % 稳定点用实心标记
            scatter3(point(1), point(2), point(3), 200, 'g', 'filled');
        else
            % 不稳定点用空心标记
            scatter3(point(1), point(2), point(3), 200, 'r', 'LineWidth', 2);
        end
        text(point(1)+0.05, point(2)+0.05, point(3)+0.05, point_names{i}, ...
             'FontSize', 10, 'FontWeight', 'bold');
    end
    
    % 绘制立方体边框
    draw_cube();
    
    xlabel('x (监管部门)');
    ylabel('y (共享方A)');
    zlabel('z (共享方B)');
    title('平衡点稳定性分析');
    legend({'稳定点 (ESS)', '不稳定点'}, 'Location', 'best');
    grid on;
    view(45, 30);
    xlim([-0.1, 1.1]);
    ylim([-0.1, 1.1]);
    zlim([-0.1, 1.1]);
    
    % 特征值实部分布
    subplot(2, 2, 2);
    plot_eigenvalue_distribution(equilibrium_points, params, point_names);
    
    % 相流图
    subplot(2, 2, 4);
    plot_phase_flow(params);
    
    sgtitle('演化博弈稳定性分析', 'FontSize', 14, 'FontWeight', 'bold');
end

function draw_cube()
    % 绘制单位立方体边框
    edges = [
        0 0 0; 1 0 0;
        1 0 0; 1 1 0;
        1 1 0; 0 1 0;
        0 1 0; 0 0 0;
        0 0 1; 1 0 1;
        1 0 1; 1 1 1;
        1 1 1; 0 1 1;
        0 1 1; 0 0 1;
        0 0 0; 0 0 1;
        1 0 0; 1 0 1;
        1 1 0; 1 1 1;
        0 1 0; 0 1 1;
    ];
    
    for i = 1:2:size(edges, 1)
        plot3([edges(i, 1), edges(i+1, 1)], ...
              [edges(i, 2), edges(i+1, 2)], ...
              [edges(i, 3), edges(i+1, 3)], ...
              'k-', 'LineWidth', 0.5);
    end
end

function plot_eigenvalue_distribution(equilibrium_points, params, point_names)
    hold on;
    
    for i = 1:size(equilibrium_points, 1)
        point = equilibrium_points(i, :);
        J = compute_jacobian(point, params);
        eigenvalues = eig(J);
        
        % 绘制特征值
        scatter(real(eigenvalues), imag(eigenvalues), 100, 'filled');
        
        % 标注平衡点
        for j = 1:length(eigenvalues)
            text(real(eigenvalues(j)), imag(eigenvalues(j)), ...
                 sprintf(' %s', point_names{i}), 'FontSize', 8);
        end
    end
    
    % 绘制稳定性边界
    xline(0, 'r--', 'LineWidth', 2);
    
    xlabel('实部 Re(λ)');
    ylabel('虚部 Im(λ)');
    title('特征值分布');
    grid on;
    legend('特征值', '稳定性边界', 'Location', 'best');
    
    hold off;
end

function plot_phase_flow(params)
    % 绘制相流图
    [X, Y] = meshgrid(0:0.1:1, 0:0.1:1);
    Z = 0.5 * ones(size(X));  % 固定z=0.5的切面
    
    U = zeros(size(X));
    V = zeros(size(Y));
    
    for i = 1:numel(X)
        state = [X(i); Y(i); Z(i)];
        dstate = evolutionary_dynamics_jacobi(0, state, params);
        U(i) = dstate(1);
        V(i) = dstate(2);
    end
    
    % 归一化箭头长度
    magnitude = sqrt(U.^2 + V.^2);
    U = U ./ (magnitude + eps);
    V = V ./ (magnitude + eps);
    
    quiver(X, Y, U, V, 0.5, 'b');
    xlabel('x (监管部门)');
    ylabel('y (共享方A)');
    title('相流图 (z=0.5切面)');
    axis equal;
    xlim([0, 1]);
    ylim([0, 1]);
    grid on;
end

function check_stability_conditions(params)
    % 检验四个稳定性条件
    fprintf('\n\n稳定性条件检验：\n');
    fprintf('%s\n', repmat('=', 1, 60));
    
    % 提取参数
    D = params.D;
    C3 = params.C3;
    a = params.a;
    n = params.n;
    C = params.C;
    k1 = params.k1;
    k2 = params.k2;
    p1 = params.p1;
    m = params.m;
    C2 = params.C2;
    
    % 条件1
    cond1 = (D < C3) && (a*n - n - C < 0);
    fprintf('条件1: D < C3 且 an - n - C < 0\n');
    fprintf('  D = %.1f, C3 = %.1f, an - n - C = %.1f\n', D, C3, a*n - n - C);
    fprintf('  满足条件1: %s\n\n', bool2str(cond1));
    
    % 条件2
    cond2 = (n - a*n - k1*p1 < 0) && (C2 - C3 + C2*D < 0) && (n - a*n - k2*p1 < 0);
    fprintf('条件2: n - an - k1p1 < 0 且 C2 - C3 + C2*D < 0 且 n - an - k2p1 < 0\n');
    fprintf('  n - an - k1p1 = %.1f\n', n - a*n - k1*p1);
    fprintf('  C2 - C3 + C2*D = %.1f\n', C2 - C3 + C2*D);
    fprintf('  n - an - k2p1 = %.1f\n', n - a*n - k2*p1);
    fprintf('  满足条件2: %s\n\n', bool2str(cond2));
    
    % 条件3
    cond3 = (m - C - n + a*n < 0) && (C3 < D);
    fprintf('条件3: m - C - n + an < 0 且 C3 < D\n');
    fprintf('  m - C - n + an = %.1f\n', m - C - n + a*n);
    fprintf('  C3 = %.1f, D = %.1f\n', C3, D);
    fprintf('  满足条件3: %s\n\n', bool2str(cond3));
    
    % 条件4
    cond4 = (n - a*n - m - k1*p1 < 0) && (C3 - C2 - C2*D < 0) && (n - m - a*n - k2*p1 < 0);
    fprintf('条件4: n - an - m - k1p1 < 0 且 C3 - C2 - C2*D < 0 且 n - m - an - k2p1 < 0\n');
    fprintf('  n - an - m - k1p1 = %.1f\n', n - a*n - m - k1*p1);
    fprintf('  C3 - C2 - C2*D = %.1f\n', C3 - C2 - C2*D);
    fprintf('  n - m - an - k2p1 = %.1f\n', n - m - a*n - k2*p1);
    fprintf('  满足条件4: %s\n', bool2str(cond4));
end

function dstate_dt = evolutionary_dynamics_jacobi(~, state, params)
    x = state(1);
    y = state(2);
    z = state(3);
    
    % 提取参数
    k1 = params.k1;
    k2 = params.k2;
    C2 = params.C2;
    C3 = params.C3;
    D = params.D;
    p1 = params.p1;
    C = params.C;
    m = params.m;
    a = params.a;
    n = params.n;
    
    % 复制动态方程
    dx_dt = x * (1 - x) * ((C2 * D * y * z) + D - D * y * z - C3 + C2 * z * y);
    dy_dt = y * (1 - y) * (k1 * p1 * z + C * z + m * x + a * n - n - C);
    dz_dt = z * (1 - z) * (a * n + C * y + m * x + k2 * p1 * y - n - C);
    
    dstate_dt = [dx_dt; dy_dt; dz_dt];
end

function str = bool2str(val)
    if val
        str = '是';
    else
        str = '否';
    end
end