%%
%2024.8.22

%模型以长方体几何中心建立坐标系，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_ver8_20240822()
% 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
gamma_T=0.12;% Tx速度方向的方位角 Azimuth angle of Tx speed direction 
gamma_R=0;% Rx速度方向的方位角 Azimuth angle of Rx speed direction
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=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;%隧道壁上有效散射体数量

delta_T=0; %Tx天线元素之间间距 %2024.8.22 需要把天线元素间距置0来获取ACF
delta_R=0; %Rx天线元素之间间距
%%
%仿真参数
% 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

%确定发射器和接收器车辆的初始位置,确定天线元素p1q1p2q2的位置
% 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);

%预分配时空相关函数的矩阵以加快计算速度
% Preallocate STCF matrix
STCF_LoS_needtosum = zeros(num_samples,N_effectivescatterers);
STCF_SB_1_needtointegral_needtosum = zeros(num_samples,N_effectivescatterers);
STCF_SB_1_needtosum = zeros(num_samples,N_effectivescatterers);
STCF_SB_2_needtointegral_needtosum = zeros(num_samples,N_effectivescatterers);
STCF_SB_2_needtosum = zeros(num_samples,N_effectivescatterers);
STCF_SB_3_needtointegral_needtosum = zeros(num_samples,N_effectivescatterers);
STCF_SB_3_needtosum = zeros(num_samples,N_effectivescatterers);



%计算每个时间样本的CIR
% Calculate CIR for each time sample
tau = 0 : 10^(-2)/(num_samples-1) : 10^(-2);%时延成为自变量

t=2;% 2024.8.22 由于绘的是STCF关于时空相关函数时延tau的图因此t为常量，在特定的t下变化，t不再是变量，取值为0,0.05,0.1(考虑到隧道长度10m，Tx位置，Tx速度20m/s，t最多是sqrt(2)/8=0.177s)
for index = 1:num_samples
    %重申坐标系的方向，正方向x朝右，y朝自己，z朝上，p1p2，q1q2的位置看论文绘图，p1相对于Tx的x坐标为负，y坐标为正，z坐标为负，p2反之
p_1_initial_position = Tx_initial_position + [-1/2*delta_T*cos(psi_T)*cos(theta_T), 1/2*delta_T*cos(psi_T)*sin(theta_T), -1/2*delta_T*sin(psi_T)];
p_2_initial_position = Tx_initial_position + [1/2*delta_T*cos(psi_T)*cos(theta_T), -1/2*delta_T*cos(psi_T)*sin(theta_T), 1/2*delta_T*sin(psi_T)];
q_1_initial_position = Rx_initial_position + [-1/2*delta_R*cos(psi_R)*cos(theta_R), 1/2*delta_T*cos(psi_R)*sin(theta_R), -1/2*delta_T*sin(psi_R)];
q_2_initial_position = Rx_initial_position + [1/2*delta_R*cos(psi_R)*cos(theta_R), -1/2*delta_T*cos(psi_R)*sin(theta_R), 1/2*delta_T*sin(psi_R)];

p_2_initial_position = p_1_initial_position;%2024.8.22 需要把p2=p1，q2=q1，那么STCF就是算的就是E｛h_p_1q_1(t)*h_p_1q_1(t+tau)｝为ACF自相关函数
q_2_initial_position = q_1_initial_position;

    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
    
    %更新发射器和接收器车辆的位置 %更新各天线元素的位置   gamma是与x轴正方向所成的夹角，如果朝右上角开gamma为正，如果朝右下角开gamma为负 
    % Update positions of the transmitter and receiver vehicles
    Tx_position = Tx_initial_position + [v_T * cos(gamma_T) * t, - v_T * sin(gamma_T) * t, 0];%2024.8.21 汽车不再是直线行驶而是呈45度夹角
    Rx_position = Rx_initial_position + [v_R * cos(gamma_R) * t, - v_R * sin(gamma_R) * t, 0];
    
    p_1_position = p_1_initial_position + [v_T * cos(gamma_T) * t, - v_T * sin(gamma_T) * t, 0];
    p_2_position = p_2_initial_position + [v_T * cos(gamma_T) * t, - v_T * sin(gamma_T) * t, 0];
    q_1_position = q_1_initial_position + [v_R * cos(gamma_R) * t, - v_R * sin(gamma_R) * t, 0];
    q_2_position = q_2_initial_position + [v_R * cos(gamma_R) * t, - v_R * sin(gamma_R) * t, 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);
        
        %计算天线元素p_1,p_2,q_1,q_2和各个散射体的距离，用于STCF的计算
        distance_p_1_S_Tx = sqrt((p_1_position(1) - scatterers_positions(1, 1))^2 + (p_1_position(2) - scatterers_positions(1, 2))^2 + (p_1_position(3) - scatterers_positions(1, 3))^2);
        distance_p_1_S_Rx = sqrt((p_1_position(1) - scatterers_positions(2, 1))^2 + (p_1_position(2) - scatterers_positions(2, 2))^2 + (p_1_position(3) - scatterers_positions(2, 3))^2);
        distance_p_1_S_tunnel= sqrt((p_1_position(1) - scatterers_positions(3, 1))^2 + (p_1_position(2) - scatterers_positions(3, 2))^2 + (p_1_position(3) - scatterers_positions(3, 3))^2);
        
        distance_p_2_S_Tx = sqrt((p_2_position(1) - scatterers_positions(1, 1))^2 + (p_2_position(2) - scatterers_positions(1, 2))^2 + (p_2_position(3) - scatterers_positions(1, 3))^2);
        distance_p_2_S_Rx = sqrt((p_2_position(1) - scatterers_positions(2, 1))^2 + (p_2_position(2) - scatterers_positions(2, 2))^2 + (p_2_position(3) - scatterers_positions(2, 3))^2);
        distance_p_2_S_tunnel= sqrt((p_2_position(1) - scatterers_positions(3, 1))^2 + (p_2_position(2) - scatterers_positions(3, 2))^2 + (p_2_position(3) - scatterers_positions(3, 3))^2);
        
        distance_q_1_S_Tx = sqrt((q_1_position(1) - scatterers_positions(1, 1))^2 + (q_1_position(2) - scatterers_positions(1, 2))^2 + (q_1_position(3) - scatterers_positions(1, 3))^2);
        distance_q_1_S_Rx = sqrt((q_1_position(1) - scatterers_positions(2, 1))^2 + (q_1_position(2) - scatterers_positions(2, 2))^2 + (q_1_position(3) - scatterers_positions(2, 3))^2);
        distance_q_1_S_tunnel= sqrt((q_1_position(1) - scatterers_positions(3, 1))^2 + (q_1_position(2) - scatterers_positions(3, 2))^2 + (q_1_position(3) - scatterers_positions(3, 3))^2);
        
        distance_q_2_S_Tx = sqrt((q_2_position(1) - scatterers_positions(1, 1))^2 + (q_2_position(2) - scatterers_positions(1, 2))^2 + (q_2_position(3) - scatterers_positions(1, 3))^2);
        distance_q_2_S_Rx = sqrt((q_2_position(1) - scatterers_positions(2, 1))^2 + (q_2_position(2) - scatterers_positions(2, 2))^2 + (q_2_position(3) - scatterers_positions(2, 3))^2);
        distance_q_2_S_tunnel= sqrt((q_2_position(1) - scatterers_positions(3, 1))^2 + (q_2_position(2) - scatterers_positions(3, 2))^2 + (q_2_position(3) - scatterers_positions(3, 3))^2);
        

        %计算公式中的夹角(方位角)以上标为向量终点，下标为向量起点，天线以p2，q2为向量终点，p1，q1为向量起点
        %Calculate the angles in the formulas
        
        %LoS路径中的alpha_LoS_T, beta_LoS_T, alpha_LoS_R,beta_LoS_R， 由于车子都是直线同向行驶可简化
        alpha_LoS_T=theta_T;
        beta_LoS_T=psi_T;
        alpha_LoS_R=theta_R;
        beta_LoS_R=psi_R;
%%
%算Tx和Tx处散射体，Rx处散射体和隧道壁散射体所成方位角和仰角   
        %SB_1路径中的角度，端点分别为Tx，Rx中央那点以及各个散射体点，求方位角时用投影向量和x轴单位向量做点积，求仰角时用向量和z轴单位向量做点积，注意正负号和是theta还是pi-theta
        vector_n_1_T = scatterers_positions(1,:) - Tx_position;%Tx处散射体和Tx所成向量，散射体为终点，Tx为起点
        project_n_1_T_xyplane = [ vector_n_1_T(1), vector_n_1_T(2), 0 ];% 散射体n_1-Tx连线向量在x-y平面上的投影
        %对于方位角来说，范围是[-pi,pi]，所以无需其他角度转换
        alpha_n_1_T = acos((project_n_1_T_xyplane*[1,0,0]')/(norm(project_n_1_T_xyplane)*1));
        
        %对于仰角来讲，范围是[-pi/2,pi/2]，所以在向量的箭头朝着z轴负平面时得进行角度转换，若朝着z轴正平面，经过点积计算出来的是pi/2-beta，所以只需要用pi/2-theta即可得到pi/2-(pi/2-beta)=beta
        if (vector_n_1_T(3)) > 0 %判断向量是在z轴上方还是下方
           beta_n_1_T = pi/2 - acos((vector_n_1_T*[0,0,1]')/(norm(vector_n_1_T)*1));
        else
           beta_n_1_T = acos((vector_n_1_T*[0,0,1]')/(norm(vector_n_1_T)*1)) - pi/2;
        end
        
        
        
        vector_n_2_T = scatterers_positions(2,:) - Tx_position;%算n_2散射体要用scatterers_positions的第二行数据
        project_n_2_T_xyplane = [ vector_n_2_T(1), vector_n_2_T(2), 0 ];
        alpha_n_2_T = acos((project_n_2_T_xyplane*[1,0,0]')/(norm(project_n_2_T_xyplane)*1));       
        if (vector_n_2_T(3)) > 0
           beta_n_2_T = pi/2 - acos((vector_n_2_T*[0,0,1]')/(norm(vector_n_2_T)*1));
        else
           beta_n_2_T = acos((vector_n_2_T*[0,0,1]')/(norm(vector_n_2_T)*1)) - pi/2;
        end

        
        
        vector_n_3_T = scatterers_positions(3,:) - Tx_position;%算n_3散射体要用scatterers_positions的第三行数据
        project_n_3_T_xyplane = [ vector_n_3_T(1), vector_n_3_T(2), 0 ];
        alpha_n_3_T = acos((project_n_3_T_xyplane*[1,0,0]')/(norm(project_n_3_T_xyplane)*1));       
        if (vector_n_3_T(3)) > 0
           beta_n_3_T = pi/2 - acos((vector_n_3_T*[0,0,1]')/(norm(vector_n_3_T)*1));
        else
           beta_n_3_T = acos((vector_n_3_T*[0,0,1]')/(norm(vector_n_3_T)*1)) - pi/2;
        end

 %开始算Rx和Tx处散射体，Rx处散射体和隧道壁散射体所成方位角和仰角    
        vector_n_1_R = scatterers_positions(1,:) - Rx_position;
        project_n_1_R_xyplane = [ vector_n_1_R(1), vector_n_1_R(2), 0 ];
       
        alpha_n_1_R = acos((project_n_1_R_xyplane*[1,0,0]')/(norm(project_n_1_R_xyplane)*1));
        
       
        if (vector_n_1_R(3)) > 0 
           beta_n_1_R = pi/2 - acos((vector_n_1_R*[0,0,1]')/(norm(vector_n_1_R)*1));
        else
           beta_n_1_R = acos((vector_n_1_R*[0,0,1]')/(norm(vector_n_1_R)*1)) - pi/2;
        end
        
        
        
        vector_n_2_R = scatterers_positions(2,:) - Rx_position;
        project_n_2_R_xyplane = [ vector_n_2_R(1), vector_n_2_R(2), 0 ];
        alpha_n_2_R = acos((project_n_2_R_xyplane*[1,0,0]')/(norm(project_n_2_R_xyplane)*1));       
        if (vector_n_2_R(3)) > 0
           beta_n_2_R = pi/2 - acos((vector_n_2_R*[0,0,1]')/(norm(vector_n_2_R)*1));
        else
           beta_n_2_R = acos((vector_n_2_R*[0,0,1]')/(norm(vector_n_2_R)*1)) - pi/2;
        end

        
        
        vector_n_3_R = scatterers_positions(3,:) - Rx_position;
        project_n_3_R_xyplane = [ vector_n_3_R(1), vector_n_3_R(2), 0 ];
        alpha_n_3_R = acos((project_n_3_R_xyplane*[1,0,0]')/(norm(project_n_3_R_xyplane)*1));       
        if (vector_n_3_R(3)) > 0
           beta_n_3_R = pi/2 - acos((vector_n_3_R*[0,0,1]')/(norm(vector_n_3_R)*1));
        else
           beta_n_3_R = acos((vector_n_3_R*[0,0,1]')/(norm(vector_n_3_R)*1)) - pi/2;
        end
  
        
        
%%        
        % 计算没有经过有效散射体加总的CIR，无穷加总做不到只能进行有限求和
        % Calculating CIR without effective scatterer summation, infinite summation cannot be achieved, only finite summation can be performed
        CIR_LoS_needtosum(index,N_effectivescatterers) = sqrt(K/(K+1))*exp(1i*phi_0)*exp(-1i*2*pi/lambda*distance_Tx_Rx)*exp(1i*2*pi/lambda*(v_T*t*cos(alpha_LoS_T-gamma_T)*cos(beta_LoS_T)+v_R*t*cos(alpha_LoS_R-gamma_R)*cos(beta_LoS_R)));
    
        CIR_SB_1_needtosum(index,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))*exp(1i*2*pi/lambda*(v_T*t*cos(alpha_n_1_T-gamma_T)*cos(beta_n_1_T)))*exp(1i*2*pi/lambda*(v_R*t*cos(alpha_n_1_R-gamma_R)*cos(beta_n_1_R)));
        CIR_SB_2_needtosum(index,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))*exp(1i*2*pi/lambda*(v_T*t*cos(alpha_n_2_T-gamma_T)*cos(beta_n_2_T)))*exp(1i*2*pi/lambda*(v_R*t*cos(alpha_n_2_R-gamma_R)*cos(beta_n_2_R)));
        CIR_SB_3_needtosum(index,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))*exp(1i*2*pi/lambda*(v_T*t*cos(alpha_n_3_T-gamma_T)*cos(beta_n_3_T)))*exp(1i*2*pi/lambda*(v_R*t*cos(alpha_n_3_R-gamma_R)*cos(beta_n_3_R)));
    
%%
%时空相关函数仿真参数


%计算天线元素之间距离
        distance_p_1_q_1= sqrt((p_1_position(1) -q_1_position(1))^2 + (p_1_position(2) - q_1_position(2))^2 + (p_1_position(3) - q_1_position(3))^2);
        distance_p_2_q_2= sqrt((p_2_position(1) -q_2_position(1))^2 + (p_2_position(2) - q_2_position(2))^2 + (p_2_position(3) - q_2_position(3))^2);
        
%计算时空相关函数STCF_LoS
STCF_LoS_needtosum(index,N_effectivescatterers) = K/(K+1)*exp(1i*2*pi/lambda*(distance_p_1_q_1-distance_p_2_q_2))*exp(1i*2*pi/lambda*tau(index)*(v_T*cos(alpha_LoS_T-gamma_T)-v_R*cos(alpha_LoS_R-gamma_R)));

%计算时空相关函数STCF_SB
        distance_p_1_q_1_n_1 = distance_p_1_S_Tx + distance_q_1_S_Tx;%Tx处散射体到p_1, q_1两天线元素之间的距离
        distance_p_2_q_2_n_1 = distance_p_2_S_Tx + distance_q_2_S_Tx;
        
        distance_p_1_q_1_n_2 = distance_p_1_S_Rx + distance_q_1_S_Rx;%Rx处散射体到p_1, q_1两天线元素之间的距离
        distance_p_2_q_2_n_2 = distance_p_2_S_Rx + distance_q_2_S_Rx;
        distance_p_1_q_1_n_3 = distance_p_1_S_tunnel + distance_q_1_S_tunnel;%Tunnel处散射体到p_1, q_1两天线元素之间的距离
        distance_p_2_q_2_n_3 = distance_p_2_S_tunnel + distance_q_2_S_tunnel;
        
%除了积分上下限，概率密度函数和积分符号之外的STCF_SB_1公式
STCF_SB_1_needtointegral_needtosum(index,N_effectivescatterers) = eta_SB_1/(K+1)*exp(1i*2*pi/lambda*(distance_p_1_q_1_n_1-distance_p_2_q_2_n_1))*exp(1i*2*pi/lambda*v_T*tau(index)*cos(alpha_n_1_T-gamma_T)*cos(beta_n_1_T))*exp(-1i*2*pi/lambda*v_R*tau(index)*cos(alpha_n_1_R-gamma_R)*cos(beta_n_1_R));
%对于将定积分离散化的情况，这里直接使用随机生成的散射体，sigma(f(x_i)*delta(x_i))可利用散射体的冯米塞斯概率密度函数*（1/散射体的数目）来近似
STCF_SB_1_needtosum(index,N_effectivescatterers) = STCF_SB_1_needtointegral_needtosum(1,N_effectivescatterers) .* vonmisesalpha(k_1,alpha_n_1_T,alpha_mu_1) .* vonmisesbeta(k_2,beta_n_1_T,alpha_mu_2) * 1/N_effectivescatterers * 1/N_effectivescatterers; 

% %定积分离散化求和 (这一版弃了，直接用前面随机生成的散射体来离散求和，只要点够多就可以近似定积分)
% number_for_integral = 100; %用于求离散积分的取样点数
% alpha_integral= -pi : 2*pi/(number_for_integral-1) : pi; %用于求离散积分的方位角取值范围及点
% beta_integral= -pi/2 : pi/(number_for_integral-1) : pi/2; %用于求离散积分的方位角取值范围及点

%STCF_SB_2
STCF_SB_2_needtointegral_needtosum(index,N_effectivescatterers) = eta_SB_2/(K+1)*exp(1i*2*pi/lambda*(distance_p_1_q_1_n_2-distance_p_2_q_2_n_2))*exp(1i*2*pi/lambda*v_T*tau(index)*cos(alpha_n_2_T-gamma_T)*cos(beta_n_2_T))*exp(-1i*2*pi/lambda*v_R*tau(index)*cos(alpha_n_2_R-gamma_R)*cos(beta_n_2_R));

STCF_SB_2_needtosum(index,N_effectivescatterers) = STCF_SB_2_needtointegral_needtosum(1,N_effectivescatterers) .* vonmisesalpha(k_1,alpha_n_2_R,alpha_mu_1) .* vonmisesbeta(k_2,alpha_n_2_R,alpha_mu_2) * 1/N_effectivescatterers * 1/N_effectivescatterers; 

%STCF_SB_3
STCF_SB_3_needtointegral_needtosum(index,N_effectivescatterers) = eta_SB_3/(K+1)*exp(1i*2*pi/lambda*(distance_p_1_q_1_n_3-distance_p_2_q_2_n_3))*exp(1i*2*pi/lambda*v_T*tau(index)*cos(alpha_n_3_T-gamma_T)*cos(beta_n_3_T))*exp(-1i*2*pi/lambda*v_R*tau(index)*cos(alpha_n_3_R-gamma_R)*cos(beta_n_3_R));

%隧道壁上的散射体为均匀分布，应用均匀分布的概率密度分布函数
[PDF_x_tunnelwall,PDF_y_tunnelwall,PDF_z_tunnelwall] = uniformtunnelscatters(tunnel_length,tunnel_width,tunnel_height);
STCF_SB_3_needtosum(1,N_effectivescatterers) = STCF_SB_3_needtointegral_needtosum(1,N_effectivescatterers) .* PDF_x_tunnelwall .* PDF_y_tunnelwall * 1/N_effectivescatterers * 1/N_effectivescatterers; 

    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

%计算STCF的贡献总和 （STCF_LoS，STCF_SB_1, STCF_SB_2, STCF_SB_3）

STCF_LoS = sum(STCF_LoS_needtosum,2);
STCF_SB_1 = sum(STCF_SB_1_needtosum,2);
STCF_SB_2 = sum(STCF_SB_2_needtosum,2);
STCF_SB_3 = sum(STCF_SB_3_needtosum,2);
STCF_total = STCF_LoS + STCF_SB_1 + STCF_SB_2 + STCF_SB_3;



% %%
% % Plot the CIR with time t
% figure;
% plot(time,abs(CIR_total))%用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}');
% legend(num2str(N_effectivescatterers)) %显示不同的散射体数目
% 
% % title('Channel Impulse Response (CIR) with different Recian factors');
% % legend(num2str(K)) % 显示不同的莱斯因子
% 
% % title('Channel Impulse Response (CIR) with different amplitude gains');
% % legend('0,0,0') % 显示不同的信道增益
% 
% xlabel('Time (s)');
% ylabel('CIR');
% 
% 
% 
%%
%绘出STCF的图
%Plot the STCF with time t

% x=delta_T/lambda;
% y=(abs(STCF_total))';
% 
% p = polyfit(x, y, 10);
% xfit = 0:10/99:10;
% yfit = polyval(p, xfit);
% % plot(xfit, yfit, 'r-'); 
% hold on
plot(tau/10^(-3),abs(STCF_total),'-g*')
title('STCF_{total} with different moving time t')
xlabel('$$\tau(ms)$$','Interpreter','latex');%标签打复杂数学公式需要按此格式
ylabel('STCF')
hold on
% %%
% %绘出多普勒频谱
% PSD = fft(STCF_total);
% % PSD_symmetry = flipud(PSD);
% plot(1./time, 20*log10(abs(PSD)))
% title('PSD')
% hold on;




% %%
% % 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_i)%计算方位角服从冯米塞斯分布的函数 k,alpha_mu为常数,alpha_i为变量 besseli(0,k_i)为以k_i为变量的0阶贝塞尔修正函数
vonmisesPDF=exp(k_i.*cos(alpha_i-alpha_mu_i))./(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_i)%计算仰角服从冯米塞斯分布的函数 k,beta_mu为常数,beta_i为变量 besseli(0,k_i)为以k_i为变量的0阶贝塞尔修正函数
vonmisesPDF=exp(k_i.*cos(beta_i-beta_mu_i))./(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

