%% 参数设置
clear; clc;
fc = 3.5e9;                  % 载波频率 28 GHz
c = physconst('LightSpeed');
lambda = c/fc;
d = lambda/2;               % 天线间距
N = 32;                    % 基站天线数量
v = 20;                     % 用户速度 (m/s)
delta_t = 0.625e-3;         % 时间间隔 (s)
num_time_points = 7;        % 时间点数
max_reflections = 2;        % 最大反射次数

%% 用户初始位置生成
rng('shuffle');
valid_position = false;
while ~valid_position
    x0 = 50*(rand() - 0.5); % x方向 [-25m, 25m]
    z0 = 50*rand() + 10;    % z方向 [10m, 60m]
    if sqrt(x0^2 + z0^2) > 15
        valid_position = true;
    end
end

%% 创建反射面模型
% 定义两个反射面参数（东侧和西侧）
reflector_positions = [20, -15]; % 反射面x坐标
reflector_rcs = [10, 5];         % 雷达截面积 (m²)

%% 创建收发设备
ula = phased.ULA('NumElements', N, 'ElementSpacing', d);
tx_pos = [0; 0; 0]; % 基站位置

%% 主循环
H = zeros(num_time_points, N);

for t_idx = 1:num_time_points
    % 更新用户位置
    current_z = z0 + v*(t_idx-1)*delta_t;
    rx_pos = [x0; 0; current_z];
    
    % 生成多径参数（LOS + 反射路径）
    [paths, angles, delays] = generate_multipath(tx_pos, rx_pos, reflector_positions, reflector_rcs, lambda, c);
    
    % 生成信道响应
    arrayResponse = zeros(N, 1);
    for p = 1:length(paths)
        % 路径参数
        theta = angles(p);
        alpha = paths(p).gain * exp(-1j*2*pi*fc*delays(p));
        
        % 多普勒计算
        dopplerShift = (v/lambda)*cosd(theta);
        phaseShift = exp(1j*2*pi*dopplerShift*(t_idx-1)*delta_t);
        
        % 阵列响应
        steeringVector = exp(-1j*2*pi*d*(0:N-1)'*sind(theta)/lambda);
        arrayResponse = arrayResponse + alpha * phaseShift * steeringVector;
    end
    
    H(t_idx, :) = arrayResponse.';
end
save('H_raytracing.mat','H');
% 可视化
figure;
subplot(2,1,1);
plot(1:N, abs(H), 'LineWidth', 1.5);
title('Channel Magnitude Response');
xlabel('Antenna Index'); ylabel('Magnitude');

subplot(2,1,2);
plot(1:N, unwrap(angle(H)), 'LineWidth', 1.5);
title('Channel Phase Distribution');
xlabel('Antenna Index'); ylabel('Phase (rad)');
%% 辅助函数：生成多径参数（修正版）
function [paths, angles, delays] = generate_multipath(tx_pos, rx_pos, reflector_x, rcs_values, lambda, c)
    % LOS路径计算
    los_vec = rx_pos - tx_pos;
    los_dist = norm(los_vec);
    paths(1).gain = 1/(4*pi*los_dist/lambda); % 自由空间路径损耗
    angles(1) = atan2d(los_vec(1), los_vec(3)); % 方位角
    delays(1) = los_dist/c;
    
    % 反射路径计算
    for r = 1:length(reflector_x)
        ref_x = reflector_x(r);
        % 镜像法计算反射路径
        image_pos = [2*ref_x - tx_pos(1); tx_pos(2); tx_pos(3)];
        ref_vec = rx_pos - image_pos;
        ref_dist = norm(ref_vec);
        
        % 反射增益计算（含RCS影响）
        reflection_coeff = sqrt(rcs_values(r)/(4*pi)) / (4*pi*ref_dist/lambda);
        paths(r+1).gain = reflection_coeff;
        angles(r+1) = atan2d(ref_vec(1), ref_vec(3));
        delays(r+1) = ref_dist/c;
    end
end

