%%
%模型以长方体几何中心建立坐标系，Tx初始位置位于长方体隧道的下车道最左侧中央（总共两个车道），Tx坐标为(x,y,z)=[0,1/4*tunnel_width,H_T/2],Rx坐标为(x,y,z)=[D,1/4*tunnel_width,H_R/2]
%The model establishes a coordinate system based on the geometric center of a rectangular prism. The initial position of Tx is located at the center of the leftmost lane of the rectangular tunnel (two lanes in total). The Tx coordinates are (x, y, z)=[0,1/4 * tunnelwidth, H_T/2], and the Rx coordinates are (x, y, z)=[D, 1/4 * tunnelwidth, H_R/2]
%两个车辆行驶方向一致，并且都是直线行驶，然后天线的角度，高度全都一样，写信道冲激响应那些公式就很简单不用代那么多参数
%Two vehicles are traveling in the same direction and in a straight line, and the angle and height of the antenna are all the same. Writing the formulas for channel impulse response is simple and does not require so many parameters
function Master_Project_UWTSD_matlab_code_ver4_20240816()
% Define the tunnel and vehicle parameters
tunnel_length = 50; % 隧道长度，单位是米 Tunnel length, in meters
tunnel_width = 10; % 隧道宽度
tunnel_height = 5; % 隧道高度
num_scatterers=3;%散射体分布有三处，Tx处，Rx处和隧道壁上，1代表Tx圆柱壁上，2代表Rx圆柱壁上，3代表隧道壁上 There are three locations where scatterers are distributed: Tx, Rx, and on the tunnel wall. 1 represents on the Tx cylindrical wall, 2 represents on the Rx cylindrical wall, and 3 represents on the tunnel wall

phi_0=0;%信号路径计算时的初始相位 Initial phase during signal path calculation

v_T = 20; % Tx的速度 speed of transmitter vehicle in m/s
v_R = 15; % Rx的速度 speed of receiver vehicle in m/s
f_c = 5.9e9; % 汽车通信载频 5.9 GHz for V2V communication
lambda = 3e8 / f_c;% 波长 wavelength
R_T= 3;% Tx处圆柱半径 Cylinder radius at Tx
R_R= 3;% Rx处圆柱半径 Cylinder radius at Rx
H_T= 1;% Tx处圆柱高度 Cylinder height at Tx
H_R= 1;% Rx处圆柱高度 Cylinder height at Rx
D= 25; % Tx和Rx之间的距离 Distance between Tx and Rx
theta_T=pi/2;% Tx天线的方位角 Azimuth angle of Tx antenna
theta_R=pi/2;% Rx天线的方位角 Azimuth angle of Rx antenna
psi_T=pi/4;% Tx天线的仰角 Elevation angle of Tx antenna
psi_R=pi/4;% Rx天线的仰角 Elevation angle of Rx antenna
alpha_mu_1=pi/4;% 散射体方位角服从的冯米塞斯分布方位角平均值 The mean azimuth value of the von Mises distribution followed by the scattering body's azimuth angle
alpha_mu_2=0;% 散射体仰角服从的冯米塞斯分布仰角平均值 The average elevation angle of the von Mises distribution followed by the scatterer's elevation angle
K=0.5;% 莱斯因子K Ricean factor K
k_1=5;% 散射体方位角服从的冯米塞斯分布的集中度参数k_1 The concentration parameter k_1 of the von Mises distribution followed by the azimuth angle of the scatterer
k_2=1;% 散射体仰角服从的冯米塞斯分布的集中度参数k_2 The concentration parameter k_2 of the von Mises distribution followed by the elevation angle of the scatterer
eta_SB_1=0.3;% SB路径Tx-S_T-Rx的信道振幅增益 Channel amplitude gain of SB path Tx-S_T-Rx
eta_SB_2=0.3;% SB路径Tx-S_R-Rx的信道振幅增益 Channel amplitude gain of SB path Tx-S_R-Rx
eta_SB_3=0.4;% SB路径Tx-S_Tunnel-Rx的信道振幅增益 Channel amplitude gain of SB path Tx-S_Tunnel-Rx

N_effectivescatterers=10000; %设置有效散射体数量统一以便for循环数一样
N_Tx=100;%Tx有效散射体数量
N_Rx=100;%Rx有效散射体数量
N_tunnel=100;%隧道壁上有效散射体数量

%%
%仿真参数
% Simulation parameters
num_samples = 100; % 采样点数 Number of time samples
sampling_interval = 0.001; % 采样间距 Time interval between samples (s)
time = (0:num_samples-1) * sampling_interval; %采样时刻 sample time

%确定发射器和接收器车辆的初始位置
% Define the initial positions of the transmitter and receiver vehicles
tx_initial_position = [0, tunnel_width / 4, tunnel_height / 2]; % 位置矩阵维度为1*3，分别代表x，y，z。 Transmitter starts at the beginning of the tunnel 
rx_initial_position = [D, tunnel_width / 4, tunnel_height / 2]; % Tx，Rx距离为D。 Receiver starts D meters from Tx

%%
%预分配信道冲激响应的矩阵以加快计算速度
% Preallocate CIR matrix
CIR_LoS_needtosum = zeros(num_samples,N_effectivescatterers);%此处分配的是需要求和的CIR因为公式中要求对同一时间t的所有散射体信道冲激响应求和 The CIR assigned here requires summation because the formula requires summation of the impulse responses of all scatterer channels at the same time t
CIR_SB_1_needtosum = zeros(num_samples,N_effectivescatterers);
CIR_SB_2_needtosum = zeros(num_samples,N_effectivescatterers);
CIR_SB_3_needtosum = zeros(num_samples,N_effectivescatterers);

%计算每个时间样本的CIR
% Calculate CIR for each time sample
for t = 1:num_samples
    for n_i = 1:N_effectivescatterers %对同一时间t的不同散射体下的CIR进行加总，公式上写的是无穷但是仿真只能做有限的 The CIR under different scatterers at the same time t is summed up, and the formula is infinite, but simulation can only do finite calculations
    
    %更新发射器和接收器车辆的位置    
    % Update positions of the transmitter and receiver vehicles
    Tx_position = tx_initial_position + [v_T * time(t), 0, 0];
    Rx_position = rx_initial_position + [v_R * time(t), 0, 0];
    
    %随机生成3个不同区域隧道内散射体的位置
    % Randomly generate the positions of scatterers within the tunnel in 3 different areas
    %隧道壁上的散射体不考虑天花板和地面的，只考虑左右两侧的，简化运算这样我可以固定纵坐标成两个数,-1/2*tunnel_width和1/2*tunnel_width，randi函数用来取这两个中的一个随机值
    %The scatterers on the tunnel wall do not consider the ceiling and ground, only the left and right sides. Simplify the calculation so that I can fix the vertical axis into two numbers, -1/2 * tunnelwidth and 1/2 * tunnelwidth. The randi function is used to take a random value from one of these two
    scatterers_positions_tunnelwall_y=[-1/2 * tunnel_width, 1/2 * tunnel_width];
    %后续对于散射体分布位置的优化可以改成分布于Tx和Rx之间，排除掉一些不符合几何单射路径的情况，但是计算会比较复杂
    %The optimization of the distribution position of scatterers in the future can be changed to be distributed between Tx and Rx, eliminating some situations that do not conform to the geometric single shot path, but the calculation will be more complex
    scatterers_positions = [Tx_position(1) + R_T * (2*rand()-1), Tx_position(2) + R_T * (2*rand()-1), Tx_position(3) + H_T * (2*rand()-1);% Tx处散射体坐标 Coordinate of scatterer at Tx location
                            Rx_position(1) + R_R * (2*rand()-1), Rx_position(2) + R_R * (2*rand()-1), Rx_position(3) + H_R * (2*rand()-1);% Rx处散射体坐标 Coordinate of scatterer at Rx location
                            tunnel_length * rand(),scatterers_positions_tunnelwall_y(randi(2)),tunnel_height * rand()];% 隧道壁上散射体坐标 Coordinates of scatterers on the tunnel wall
%%
        % 计算从Tx和Rx到散射体的距离
        % Calculate distances from Tx and Rx to the scatterers
        distance_Tx_Rx= sqrt((Tx_position(1) -Rx_position(1))^2 + (Tx_position(2) - Rx_position(2))^2 + (Tx_position(3) - Rx_position(3))^2);
        
        distance_Tx_S_Tx = sqrt((Tx_position(1) - scatterers_positions(1, 1))^2 + (Tx_position(2) - scatterers_positions(1, 2))^2 + (Tx_position(3) - scatterers_positions(1, 3))^2);
        distance_Tx_S_Rx = sqrt((Tx_position(1) - scatterers_positions(2, 1))^2 + (Tx_position(2) - scatterers_positions(2, 2))^2 + (Tx_position(3) - scatterers_positions(2, 3))^2);
        distance_Tx_S_tunnel= sqrt((Tx_position(1) - scatterers_positions(3, 1))^2 + (Tx_position(2) - scatterers_positions(3, 2))^2 + (Tx_position(3) - scatterers_positions(3, 3))^2);
        
        distance_S_Tx_Rx = sqrt((Tx_position(1) - scatterers_positions(1, 1))^2 + (Tx_position(2) - scatterers_positions(1, 2))^2 + (Tx_position(3) - scatterers_positions(1, 3))^2);
        distance_S_Rx_Rx = sqrt((Tx_position(1) - scatterers_positions(2, 1))^2 + (Tx_position(2) - scatterers_positions(2, 2))^2 + (Tx_position(3) - scatterers_positions(2, 3))^2);
        distance_S_tunnel_Rx = sqrt((Tx_position(1) - scatterers_positions(3, 1))^2 + (Tx_position(2) - scatterers_positions(3, 2))^2 + (Tx_position(3) - scatterers_positions(3, 3))^2);

        
        % 计算没有经过有效散射体加总的CIR，无穷加总做不到只能进行有限求和
        % Calculating CIR without effective scatterer summation, infinite summation cannot be achieved, only finite summation can be performed
        CIR_LoS_needtosum(t,N_effectivescatterers) = sqrt(K/(K+1))*exp(1i*phi_0)*exp(-1i*2*pi/lambda*distance_Tx_Rx);
    
        CIR_SB_1_needtosum(t,N_effectivescatterers) = sqrt(eta_SB_1/(K+1))*1/(sqrt(N_effectivescatterers))*exp(1i*phi_0)*exp(-1i*2*pi/lambda*(distance_Tx_S_Tx+distance_S_Tx_Rx));
        CIR_SB_2_needtosum(t,N_effectivescatterers) = sqrt(eta_SB_2/(K+1))*1/(sqrt(N_effectivescatterers))*exp(1i*phi_0)*exp(-1i*2*pi/lambda*(distance_Tx_S_Rx+distance_S_Rx_Rx));
        CIR_SB_3_needtosum(t,N_effectivescatterers) = sqrt(eta_SB_3/(K+1))*1/(sqrt(N_effectivescatterers))*exp(1i*phi_0)*exp(-1i*2*pi/lambda*(distance_Tx_S_tunnel+distance_S_tunnel_Rx));
    end

end

%%
% 计算CIR的贡献总和（LoS.SB_1，SB_2。sb3）
% Sum contributions from CIR (LoS. SB_1, SB_2. SB_3)
% CIR_LoS因为是LoS径没有散射体遮挡所以不需要对多个有效散射体的CIR进行求和，然而实际仿真过程中发现上面的for循环虽然LoS没有参与但是会给他赋值0，导致不是num_samples*1矩阵而是num_samples*N_effectivescatterers矩阵
% CIR_LoS does not need to sum the CIRs of multiple effective scatterers because it is a LoS path with no scatterer obstruction. However, in the actual simulation process, it was found that although LoS did not participate, the for loop above would assign it a value of 0, resulting in a num_Samples * 1 matrix instead of a num_Samples * N_effectors matrix
CIR_LoS = sum(CIR_LoS_needtosum, 2);
CIR_SB_1 = sum(CIR_SB_1_needtosum, 2);
CIR_SB_2 = sum(CIR_SB_2_needtosum, 2);
CIR_SB_3 = sum(CIR_SB_3_needtosum, 2);
CIR_total = CIR_LoS + CIR_SB_1 + CIR_SB_2 + CIR_SB_3;%叠加原理加总CIR。 Sum up CIR based on superposition principle

% Plot the CIR with time t
figure;
plot(time,abs(CIR_total),'-ks')%用abs返回复数的幅值以求CIR的幅值大小 Using ABS to return the amplitude of complex numbers to determine the magnitude of CIR
title('Channel Impulse Response (CIR) with different N_{effectivescatterers}');
xlabel('Time (s)');
ylabel('CIR');
legend(num2str(N_effectivescatterers))
% %%
% % Calculate the temporal ACF
% ACF = xcorr(CIR_total, 'coeff');
% 
% % Plot the ACF
% figure;
% plot((-num_samples+1:num_samples-1) * sampling_interval, abs(ACF));
% title('Temporal Autocorrelation Function (ACF)');
% xlabel('Time Lag (s)');
% ylabel('Autocorrelation');
% 
% % Define spatial positions for calculation
% spatial_positions = linspace(0, tunnel_length, num_samples);
% 
% % Preallocate CCF array
% CCF = zeros(num_samples, 1);
% 
% % Calculate the spatial CCF
% for k = 1:num_samples
%     % Calculate the integer shift value for spatial position
%     shift_val = round((spatial_positions(k) / tunnel_length) * num_samples);
%     
%     % Ensure the shift value is within the bounds
%     shift_val = mod(shift_val, num_samples);
%     
%     % Compute the CCF using circshift with integer shift value
%     CCF(k) = mean(CIR_total .* circshift(CIR_total, shift_val));
% end
% 
% % Plot the CCF
% figure;
% plot(spatial_positions, abs(CCF));
% title('Spatial Cross-Correlation Function (CCF)');
% xlabel('Spatial Position (m)');
% ylabel('Cross-Correlation');
% 
% % Calculate Doppler spectrum
% doppler_spectrum = abs(fftshift(fft(CIR_total))).^2;
% 
% % Frequency axis
% doppler_freq = linspace(-0.5/sampling_interval, 0.5/sampling_interval, num_samples);
% 
% % Plot Doppler PSD
% figure;
% plot(doppler_freq, 10*log10(doppler_spectrum));
% title('Doppler Power Spectral Density (PSD)');
% xlabel('Doppler Frequency (Hz)');
% ylabel('Power/Frequency (dB/Hz)');
% 
% % Visualization of the tunnel, transmitter, receiver, and scattering clusters
% figure;
% hold on;
% 
% % Plot the tunnel
% fill3([0 tunnel_length tunnel_length 0], [0 0 tunnel_width tunnel_width], [0 0 0 0], [0.8 0.8 0.8], 'FaceAlpha', 0.5); % Floor
% fill3([0 tunnel_length tunnel_length 0], [0 0 tunnel_width tunnel_width], [tunnel_height tunnel_height tunnel_height tunnel_height], [0.8 0.8 0.8], 'FaceAlpha', 0.5); % Ceiling
% fill3([0 0 0 0], [0 0 tunnel_width tunnel_width], [0 tunnel_height tunnel_height 0], [0.8 0.8 0.8], 'FaceAlpha', 0.5); % Left wall
% fill3([tunnel_length tunnel_length tunnel_length tunnel_length], [0 0 tunnel_width tunnel_width], [0 tunnel_height tunnel_height 0], [0.8 0.8 0.8], 'FaceAlpha', 0.5); % Right wall
% 
% % Plot the initial positions of the transmitter and receiver vehicles
% plot3(tx_initial_position(1), tx_initial_position(2), tx_initial_position(3), 'rs', 'MarkerSize', 10, 'MarkerFaceColor', 'r', 'DisplayName', 'Transmitter (Tx)');
% plot3(rx_initial_position(1), rx_initial_position(2), rx_initial_position(3), 'gs', 'MarkerSize', 10, 'MarkerFaceColor', 'g', 'DisplayName', 'Receiver (Rx)');
% 
% % Plot the scattering clusters as blue circles
% plot3(scatterers_positions(:, 1), scatterers_positions(:, 2), scatterers_positions(:, 3), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b', 'DisplayName', 'Scattering Clusters');
% 
% % Annotate the plot
% title('Tunnel Model with Transmitter, Receiver, and Scattering Clusters');
% xlabel('Tunnel Length (m)');
% ylabel('Tunnel Width (m)');
% zlabel('Tunnel Height (m)');
% legend('Tunnel', 'Transmitter (Tx)', 'Receiver (Rx)', 'Scattering Clusters');
% grid on;
% view(3);
% axis equal;
% hold off;


end

%%
function vonmisesPDF=vonmisesalpha(k_i,alpha_i,alpha_mu)%计算方位角服从冯米塞斯分布的函数 k,alpha_mu为常数,alpha_i为变量 besseli(0,k_i)为以k_i为变量的0阶贝塞尔修正函数
vonmisesPDF=exp(k_i.*cos(alpha_i-alpha_mu))./(2.*pi.*besseli(0,k_i));
% besseli Modified Bessel function of the first kind.
%     I = besseli(NU,Z) is the modified Bessel function of the first kind,
%     I_nu(Z).  The order NU need not be an integer, but must be real.
%     The argument Z can be complex.  The result is real where Z is positive.
end

%%
function vonmisesPDF=vonmisesbeta(k_i,beta_i,beta_mu)%计算仰角服从冯米塞斯分布的函数 k,beta_mu为常数,beta_i为变量 besseli(0,k_i)为以k_i为变量的0阶贝塞尔修正函数
vonmisesPDF=exp(k_i.*cos(beta_i-beta_mu))./(2.*pi.*besseli(0,k_i));
% besseli Modified Bessel function of the first kind.
%     I = besseli(NU,Z) is the modified Bessel function of the first kind,
%     I_nu(Z).  The order NU need not be an integer, but must be real.
%     The argument Z can be complex.  The result is real where Z is positive.
end

%%
%需要考虑到一些不可能的几何反射，放在优化考虑，在一个时间t时不可能整个隧道的隧道壁都有散射体可以完成单射路径
function [f_x, f_y, f_z] =uniformtunnelscatters(lengthoftunnel,widthoftunnel,heightoftunnel)%隧道散射体一般来说符合均匀分布，通过输入隧道的长宽高来得到分布函数表达式
f_x=1/lengthoftunnel;
f_y=1/widthoftunnel;
f_z=1/heightoftunnel;
end

