% 分析信誉参数k1对演化的影响
function analyze_reputation_parameter()
    
    fprintf('\n开始信誉参数k1的敏感性分析...\n');
    
    % 基础参数设置
    base_params = struct();
    base_params.C2 = 50;
    base_params.D = 60;
    base_params.C3 = 50;
    base_params.p1 = 200;
    base_params.C = 40;
    base_params.m = 50;
    base_params.a = 0.6;
    base_params.n = 80;
    base_params.k2 = 1;
    
    % 三组k1值对应不同的信誉值Rep
    k1_groups = {
        [0.8, 1, 1.2],    % Rep = [0.4, 0.5, 0.6]
        [1, 1.2, 1.6],    % Rep = [0.5, 0.6, 0.8]
        [0.4, 0.6, 0.8]   % Rep = [0.2, 0.3, 0.4]
    };
    
    % 创建图形窗口
    figure('Position', [100, 100, 1500, 500]);
    
    for group_idx = 1:3
        subplot(1, 3, group_idx);
        hold on;
        
        k1_values = k1_groups{group_idx};
        colors = {'r-', 'g-', 'b-'};
        legend_labels = {};
        
        for k_idx = 1:3
            base_params.k1 = k1_values(k_idx);
            rep_value = k1_values(k_idx) / 2;
            
            % 运行仿真
            tspan = [0, 50];
            initial_state = [0.5; 0.5; 0.5];
            
            [t, trajectory] = ode45(@(t, state) evolutionary_dynamics_k1(t, state, base_params), ...
                                    tspan, initial_state);
            
            % 绘制3D轨迹
            plot3(trajectory(:, 1), trajectory(:, 2), trajectory(:, 3), ...
                  colors{k_idx}, 'LineWidth', 2);
            
            legend_labels{k_idx} = sprintf('k_1=%.1f (Rep=%.1f)', k1_values(k_idx), rep_value);
        end
        
        % 设置图形属性
        xlabel('x');
        ylabel('y');
        zlabel('z');
        title(sprintf('第%d组: k_1的影响', group_idx));
        legend(legend_labels, 'Location', 'best');
        grid on;
        view(45, 30);
        xlim([0, 1]);
        ylim([0, 1]);
        zlim([0, 1]);
        
        hold off;
    end
    
    sgtitle('信誉参数k_1的敏感性分析（图7-8）', 'FontSize', 14, 'FontWeight', 'bold');
    
    % 生成2D投影图
    generate_k1_2D_projections(base_params, k1_groups);
end

function generate_k1_2D_projections(base_params, k1_groups)
    figure('Position', [100, 100, 1500, 500]);
    
    for group_idx = 1:3
        subplot(1, 3, group_idx);
        hold on;
        
        k1_values = k1_groups{group_idx};
        colors = {[0.8, 0.2, 0.2], [0.2, 0.8, 0.2], [0.2, 0.2, 0.8]};
        markers = {'x', 'o', '+'};
        
        % 生成网格
        [X, Y] = meshgrid(0:0.1:1, 0:0.1:1);
        
        for k_idx = 1:3
            base_params.k1 = k1_values(k_idx);
            
            % 运行多次仿真
            num_sims = 20;
            tspan = [0, 50];
            
            for sim = 1:num_sims
                initial_state = [rand(); rand(); rand()];
                
                [t, trajectory] = ode45(@(t, state) evolutionary_dynamics_k1(t, state, base_params), ...
                                        tspan, initial_state);
                
                % 绘制轨迹终点
                final_x = trajectory(end, 1);
                final_y = trajectory(end, 2);
                final_z = trajectory(end, 3);
                
                scatter3(final_x, final_y, final_z, 50, ...
                        'MarkerEdgeColor', colors{k_idx}, ...
                        'Marker', markers{k_idx}, ...
                        'LineWidth', 2);
            end
        end
        
        % 设置图形属性
        xlabel('x');
        ylabel('y');
        zlabel('z');
        title(sprintf('第%d组: 演化终点分布', group_idx));
        grid on;
        view(45, 30);
        xlim([0, 1]);
        ylim([0, 1]);
        zlim([0, 1]);
        
        % 添加图例
        h1 = scatter3(nan, nan, nan, 50, 'MarkerEdgeColor', colors{1}, 'Marker', markers{1}, 'LineWidth', 2);
        h2 = scatter3(nan, nan, nan, 50, 'MarkerEdgeColor', colors{2}, 'Marker', markers{2}, 'LineWidth', 2);
        h3 = scatter3(nan, nan, nan, 50, 'MarkerEdgeColor', colors{3}, 'Marker', markers{3}, 'LineWidth', 2);
        
        rep_values = k1_values / 2;
        legend([h1, h2, h3], ...
               {sprintf('k_1=%.1f', k1_values(1)), ...
                sprintf('k_1=%.1f', k1_values(2)), ...
                sprintf('k_1=%.1f', k1_values(3))}, ...
               'Location', 'best');
        
        hold off;
    end
    
    sgtitle('信誉参数k_1的演化终点分布', 'FontSize', 14, 'FontWeight', 'bold');
end

function dstate_dt = evolutionary_dynamics_k1(~, 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