simulate_parallel.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   多核CPU并行版本
%    400个元素, 0.25° 精度耗时约1.4s
%   1600个元素, 0.25° 精度耗时约4.5s
%   1600个元素, 0.05° 精度耗时约100s
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clc;
close all;

degree_min = input('请输入仿真的角度精度:');
if degree_min > 1
    degree_min = 0.0625;
elseif degree_min < 0.05
    degree_min = 0.05;
end

tic;
%% 1. 读取天线的profile文件
% 文件名
filename = 'array_3.txt';
fprintf("仿真文件：%s\n", filename);
filename_for_title = strrep(filename, '_', '\_'); % 用 \_ 替换 _

% 打开文件
fid = fopen(filename, 'r');
if fid == -1
    error('无法打开文件 %s', filename);
end

% 读取工作频率 (f) 信息
freq_line = fgetl(fid); % 读取第一行
f = sscanf(freq_line, 'f %fGHz') * 1e9; % 提取频率并转换为 Hz
fprintf('工作频率: %.2f Hz\n', f);

% 得出波长和波数值
lambda = 3e8 / f; % 波长，公式：λ = c / f，单位 m
k = 2 * pi / lambda; % 波数，公式：k = 2π / λ

% 跳过第二行表头
header_line = fgetl(fid); % 第二行是表头

% 初始化数据存储
data = []; % 用于存储 [xn, yn, An, phi]

% 逐行读取数据
while ~feof(fid)
    line = fgetl(fid); % 读取每一行
    if ischar(line) && ~isempty(strtrim(line)) % 确保行有效
        % 使用 sscanf 提取数据
        values = sscanf(line, '%*d %f %f %f %f'); % 跳过第一列索引
        if numel(values) == 4 % 检查是否提取到4列数据
            data = [data; values']; % 动态增加数据行
        end
    end
end

% 关闭文件
fclose(fid);

% 提取数据
num = size(data, 1); % 动态确定天线元件数量
xn = data(:, 1); % x 坐标
yn = data(:, 2); % y 坐标
An = data(:, 3); % 激励幅值
phi = data(:, 4); % 激励相位 (角度制)

% 显示读取的结果
fprintf('天线元件数量: %d\n', num);
% disp('读取的数据:');
% disp(table((1:num)', xn, yn, An, phi, ...
%     'VariableNames', {'Index', 'Xn', 'Yn', 'An', 'Phi'}));

%% 2. 可视化天线阵列
% 计算坐标范围
x_min = min(xn) - 0.02; % 比最小值小一点
x_max = max(xn) + 0.02; % 比最大值大一点
y_min = min(yn) - 0.02; % 比最小值小一点
y_max = max(yn) + 0.02; % 比最大值大一点

% 获取各轴的实际范围
x_range = x_max - x_min;
y_range = y_max - y_min;
z_range = max(An) - min(An);


% 计算比例因子，使得 z 轴与 x 和 y 的范围在视觉上协调
z_factor = 0.1 * max(x_range, y_range) / z_range;

% 可视化天线阵列
figure('Units', 'pixels', 'Position', [100, 100, 1000, 700]);
scatter(xn, yn, 100, An, 'filled'); % 使用散点图，点的颜色代表激励幅值
colorbar; % 添加颜色条以表示激励幅值
colormap('jet'); % 使用 jet 颜色映射
title(['天线阵列的二维分布（', filename_for_title, '）']);
xlabel('X 坐标 (m)');
ylabel('Y 坐标 (m)');
grid on;
axis equal;

% 设置坐标范围
xlim([x_min, x_max]);
ylim([y_min, y_max]);

% 在每个点上标注相位信息
for i = 1:num
    text(xn(i), yn(i), sprintf('%.1f°', phi(i)), ...
         'VerticalAlignment', 'bottom', ...
         'HorizontalAlignment', 'right', ...
         'FontSize', 8, 'Color', 'black');
end

% 3D 天线阵列分布
figure('Units', 'pixels', 'Position', [100, 100, 600, 400]);
scatter3(xn, yn, An, 100, An, 'filled'); % 3D 散点图，颜色和大小与 An 对应
colorbar; % 添加颜色条表示激励幅值
colormap('jet'); % 使用 jet 颜色映射
title(['3D 天线阵列分布（', filename_for_title, '）']);
xlabel('X 坐标 (m)');
ylabel('Y 坐标 (m)');
zlabel('激励幅值 A_n');
grid on;
rotate3d on

% 3D 天线阵列分布
figure('Units', 'pixels', 'Position', [100, 100, 600, 400]);
scatter3(xn, yn, phi, 100, phi, 'filled'); % 3D 散点图，颜色和大小与 An 对应
colorbar; % 添加颜色条表示激励幅值
colormap('jet'); % 使用 jet 颜色映射
title(['3D 天线阵列分布（', filename_for_title, '）相位']);
xlabel('X 坐标 (m)');
ylabel('Y 坐标 (m)');
zlabel('激励幅值 A_n');
grid on;
rotate3d on

% 设置坐标范围
xlim([x_min, x_max]);
ylim([y_min, y_max]);

% 调整比例使 3D 图更美观
view(60, 60); % 设置观察角度，方位角 45°, 仰角 30°

%% 3. 辐射方向图， 主梁方向和方向性

% 定义方向角范围
point_num_of_1 = floor(1 / degree_min);
theta = linspace(0, pi, 180 * point_num_of_1); % θ 从 0° 到 180°，上半球和下半球
phi_angle = linspace(0, 2*pi, 360 * point_num_of_1); % φ_angle 从 0° 到 360° (方向角)

% 初始化辐射方向图
E = zeros(length(theta), length(phi_angle));


% 总任务数
total_steps = length(theta);

% 使用 parfor 并行化外层循环
parfor t = 1:total_steps
    % 临时变量（避免共享变量冲突）
    local_E = zeros(1, length(phi_angle));
    
    for p = 1:length(phi_angle)
        % 判断是否在下半球
        if theta(t) > pi/2
            local_E(p) = 0; % 下半球设置为零场
            continue;
        end

        sum_E = 0;
        for i = 1:num
            % 矢量点积计算：k * (r_i · u)
            phase_shift = k * (xn(i) * sin(theta(t)) * cos(phi_angle(p)) + ...
                               yn(i) * sin(theta(t)) * sin(phi_angle(p)));
            % 累加所有天线元的电场，使用元件的激励相位 phi(i)
            % sum_E = sum_E + An(i) * exp(1j * (phi(i) * pi / 180 + phase_shift));
            sum_E = sum_E + An(i) * exp(1j * (phi(i) + phase_shift));
        end
        % 计算方向图
        local_E(p) = abs(sum_E); % 取绝对值（幅度）
    end
    
    % 保存结果
    E(t, :) = local_E;
end




% 归一化方向图
E_max= max(E(:));
fprintf('最大幅度：Emax = %.2f\n', E_max);
E_norm = E / max(E(:));

% 转换为笛卡尔坐标以绘制 3D 图
[theta_grid, phi_grid] = meshgrid(theta, phi_angle);
X = E_norm' .* sin(theta_grid) .* cos(phi_grid);
Y = E_norm' .* sin(theta_grid) .* sin(phi_grid);
Z = E_norm' .* cos(theta_grid);

% 绘制 3D 辐射方向图
figure('Units', 'pixels', 'Position', [100, 100, 800, 600]); % 设置图片大小
surf(X, Y, Z, E_norm', 'EdgeColor', 'none'); % 绘制表面，颜色代表幅值
colormap('jet'); % 使用 jet 颜色映射
colorbar; % 添加颜色条表示幅值
title(['3D 天线阵列辐射方向图（', strrep(filename, '_', '\_'), '）']);
xlabel('X 坐标');
ylabel('Y 坐标');
zlabel('Z 坐标');
% axis equal; % 确保坐标比例一致
view(45, 30); % 设置观察角度
grid on;
rotate3d on;

% 主梁方向
[max_val, idx] = max(E_norm(:)); % 找到最大值
[max_t_idx, max_p_idx] = ind2sub(size(E_norm), idx); % 最大值对应的索引
main_beam_theta = theta(max_t_idx) * (180 / pi); % 转换为度
main_beam_phi = phi_angle(max_p_idx) * (180 / pi); % 转换为度
fprintf('主梁方向：θ = %.2f°，φ = %.2f°\n', main_beam_theta, main_beam_phi);

% 天线方向性
% 使用离散积分计算 ∫∫E^2 sinθ dθ dφ
d_theta = theta(2) - theta(1); % θ 间隔
d_phi = phi_angle(2) - phi_angle(1); % φ 间隔
P_radiated = sum(sum((E_norm.^2) .* sin(theta') * d_theta * d_phi)); % 辐射总功率
D = 4 * pi * (max_val^2) / P_radiated; % 方向性公式
fprintf('天线阵列方向性：D = %.2f\n', D);


%% 停止计时
% 停止计时并输出耗时
elapsed_time = toc;
disp(['计算完成，耗时: ', num2str(elapsed_time), ' 秒']);

%% 绘制俯仰角的极坐标图
% 生成 E_0phi (假设已定义 E_0phi)
E_0phi = E_norm(:, 1); % 假设 E_0phi 为第一列的方向图（θ从0°到180°）

% 扩展 theta 从 0° 到 360°
theta_extended = linspace(0, pi, 180 * point_num_of_1); % θ 从 0° 到 180°
phi_extended = linspace(0, 2*pi, 360 * point_num_of_1); % φ 从 0° 到 360°

% 使用对称性扩展 E_0phi 到完整的360°范围
E_0phi_full = [E_0phi; flipud(E_0phi)]; % 将 E_0phi 反转并连接到原始数据后

% 将 E_0phi_full 扩展为与 phi 相关的数据
E_0phi_full_extended = repmat(E_0phi_full, 1, 1)';
E_0phi_full_extended(180 * point_num_of_1 + 1:end) = 0;

% 绘制极坐标图
figure;
polarplot(phi_extended, E_0phi_full_extended, 'LineWidth', 2);

% 设置图标题和标签
title(['极坐标图: E_0phi（', strrep(filename, '_', '\_'), '）']);
rlim([0 max(E_0phi_full_extended)]); % 设置极坐标半径范围