clear all; clc;

tic

f=2.4*10^9;
lambda=3*10^8/f;
kappa=2*pi/lambda;

Z0=376.73;              %真空阻抗

K=8;                    %用户数量

%% 发射机面积 坐标
Ls_x=sqrt(0.25);                           %发射机板面宽度
Ls_y=sqrt(0.25);
delta_x=33;                         %x轴量化
delta_y=33;                         %y轴量化
s_center= [0 0 0]';                 %发射机中心

Vs=zeros(3,delta_x*delta_y);         %发射机表面坐标
for a=1:delta_x
    for b=1:delta_y
        Vs(1,a+delta_x*(b-1))=Ls_x/delta_x*a-Ls_x/2;
        Vs(2,a+delta_x*(b-1))=Ls_y/delta_y*b-Ls_y/2;
        Vs(3,a+delta_x*(b-1))=0;
    end
end

dS=Ls_x*Ls_y/delta_x/delta_y;       %发射机量化单元面积
NumdS=delta_x*delta_y;              %发射机量化单元数

%% 离散MIMO
Mx = max(floor(Ls_x*2/lambda),1);
My = max(floor(Ls_y*2/lambda),1);

M = Mx*My;                          %离散天线总数

S_m = zeros(3,M);
for mx = 1:Mx
    for my = 1:My
        S_m(1,mx + Mx*(my - 1)) = mx*lambda/2-Ls_x/2;
        S_m(2,mx + Mx*(my - 1)) = my*lambda/2-Ls_y/2;
    end
end
%% 傅里叶基展开
delta_u=9;                         %u轴基展开量化
delta_v=9;                         %v轴基展开量化
NumF=delta_u*delta_v;

f_base=zeros(NumF,NumdS);           %傅里叶基生成
for u=1:delta_u
    for v=1:delta_v
        for n=1:NumdS
            f_base(u+delta_u*(v-1),n)=1/sqrt(Ls_x*Ls_y)*exp(1j*2*pi*(u-ceil(delta_u/2))*(Vs(1,n)-Ls_x/2)/Ls_x)*exp(1j*2*pi*(v-ceil(delta_v/2))*(Vs(2,n)-Ls_y/2)/Ls_y);
        end
    end
end

%% 接收机面积 坐标
Vr=lambda^2/4/pi;
Psi_k=zeros(3,K);   %combiner初始化
for k=1:K
    Psi_k(:,k)=(2*rand(3,1)-1+2j*(rand(3,1)-0.5))/sqrt(6);
    Psi_k(:,k)=Psi_k(:,k)/norm(Psi_k(:,k));
end

%% 容量计算
P_T=10^(-4);
sigma_2=5.6*10^(-3)/kappa^2/Z0^2;      %噪声功率 V2/m2

%% 用户角度
% angle = [0:2:16]/180*pi;                %4个用户的夹角
distance = [0:1:20];
distance(1) = 0.001;

%%
length_L=length(distance);

Iteration = 10;

sum_rate_H_MIMO_Upper=zeros(length_L,Iteration);
sum_rate_PDM_1=zeros(length_L,Iteration);
sum_rate_PDM_2=zeros(length_L,Iteration);
sum_rate_PDM_3=zeros(length_L,Iteration);
sum_rate_PDM_4=zeros(length_L,Iteration);
sum_rate_MF_1=zeros(length_L,Iteration);

sum_rate_ditial = zeros(length_L,Iteration);

for a=1:Iteration
    parfor m=1:length_L

        sprintf('位置，第%d轮',m)

        DD = distance(m);

        L=2;   %用户中心的Z轴坐标
        user_position_1=zeros(3,K);       %用户坐标
        user_position_1(:,1)=[DD*cos(pi/4) DD*sin(pi/4) L];
        user_position_1(:,2)=[DD*cos(2*pi/4) DD*sin(2*pi/4) L];
        user_position_1(:,3)=[DD*cos(3*pi/4) DD*sin(3*pi/4) L];
        user_position_1(:,4)=[DD*cos(4*pi/4) DD*sin(4*pi/4) L];
        user_position_1(:,5)=[DD*cos(5*pi/4) DD*sin(5*pi/4) L];
        user_position_1(:,6)=[DD*cos(6*pi/4) DD*sin(6*pi/4) L];
        user_position_1(:,7)=[DD*cos(7*pi/4) DD*sin(7*pi/4) L];
        user_position_1(:,8)=[DD*cos(8*pi/4) DD*sin(8*pi/4) L];

        L=5;   %用户中心的Z轴坐标
        user_position_2=zeros(3,K);       %用户坐标
        user_position_2(:,1)=[DD*cos(pi/4) DD*sin(pi/4) L];
        user_position_2(:,2)=[DD*cos(2*pi/4) DD*sin(2*pi/4) L];
        user_position_2(:,3)=[DD*cos(3*pi/4) DD*sin(3*pi/4) L];
        user_position_2(:,4)=[DD*cos(4*pi/4) DD*sin(4*pi/4) L];
        user_position_2(:,5)=[DD*cos(5*pi/4) DD*sin(5*pi/4) L];
        user_position_2(:,6)=[DD*cos(6*pi/4) DD*sin(6*pi/4) L];
        user_position_2(:,7)=[DD*cos(7*pi/4) DD*sin(7*pi/4) L];
        user_position_2(:,8)=[DD*cos(8*pi/4) DD*sin(8*pi/4) L];

        L=10;   %用户中心的Z轴坐标
        user_position_3=zeros(3,K);       %用户坐标
        user_position_3(:,1)=[DD*cos(pi/4) DD*sin(pi/4) L];
        user_position_3(:,2)=[DD*cos(2*pi/4) DD*sin(2*pi/4) L];
        user_position_3(:,3)=[DD*cos(3*pi/4) DD*sin(3*pi/4) L];
        user_position_3(:,4)=[DD*cos(4*pi/4) DD*sin(4*pi/4) L];
        user_position_3(:,5)=[DD*cos(5*pi/4) DD*sin(5*pi/4) L];
        user_position_3(:,6)=[DD*cos(6*pi/4) DD*sin(6*pi/4) L];
        user_position_3(:,7)=[DD*cos(7*pi/4) DD*sin(7*pi/4) L];
        user_position_3(:,8)=[DD*cos(8*pi/4) DD*sin(8*pi/4) L];

        L=30;   %用户中心的Z轴坐标
        user_position_4=zeros(3,K);       %用户坐标
        user_position_4(:,1)=[DD*cos(pi/4) DD*sin(pi/4) L];
        user_position_4(:,2)=[DD*cos(2*pi/4) DD*sin(2*pi/4) L];
        user_position_4(:,3)=[DD*cos(3*pi/4) DD*sin(3*pi/4) L];
        user_position_4(:,4)=[DD*cos(4*pi/4) DD*sin(4*pi/4) L];
        user_position_4(:,5)=[DD*cos(5*pi/4) DD*sin(5*pi/4) L];
        user_position_4(:,6)=[DD*cos(6*pi/4) DD*sin(6*pi/4) L];
        user_position_4(:,7)=[DD*cos(7*pi/4) DD*sin(7*pi/4) L];
        user_position_4(:,8)=[DD*cos(8*pi/4) DD*sin(8*pi/4) L];

        theta_k=zeros(3,K,NumdS);
        temp=0;
        for n=1:NumdS
            for k=1:K
                theta_k(:,k,n)=(2*rand(3,1)-1+2j*(rand(3,1)-0.5))/sqrt(6);
                temp=temp+norm(theta_k(:,k,n))^2*dS;
            end
        end
        theta_k=theta_k*sqrt(P_T)/sqrt(temp);

        G_k = Freespace_G_k_generate(lambda,user_position_1,Vs,K,NumdS,delta_x,delta_y);
        [~,~,sum_rate_PDM_1(m,a)] = H_MIMO_precoding(K,NumdS,NumF,dS,G_k,Psi_k,theta_k,f_base,P_T,sigma_2,Ls_x,Ls_y);

        [~,~,sum_rate_PDM_1(m,a)] = H_MIMO_precoding(K,NumdS,NumF,dS,G_k,Psi_k,theta_k,f_base,P_T,sigma_2,Ls_x,Ls_y);
        [~,~,sum_rate_MF_1(m,a)] = MF_MIMO_precoding(K,NumdS,dS,G_k,Psi_k,theta_k,P_T,sigma_2);
       
        G_k = Freespace_G_k_generate(lambda,user_position_2,Vs,K,NumdS,delta_x,delta_y);
        [~,~,sum_rate_PDM_2(m,a)] = H_MIMO_precoding(K,NumdS,NumF,dS,G_k,Psi_k,theta_k,f_base,P_T,sigma_2,Ls_x,Ls_y);
        G_k = Freespace_G_k_generate(lambda,user_position_3,Vs,K,NumdS,delta_x,delta_y);
        [~,~,sum_rate_PDM_3(m,a)] = H_MIMO_precoding(K,NumdS,NumF,dS,G_k,Psi_k,theta_k,f_base,P_T,sigma_2,Ls_x,Ls_y);
        G_k = Freespace_G_k_generate(lambda,user_position_4,Vs,K,NumdS,delta_x,delta_y);
        [~,~,sum_rate_PDM_4(m,a)] = H_MIMO_precoding(K,NumdS,NumF,dS,G_k,Psi_k,theta_k,f_base,P_T,sigma_2,Ls_x,Ls_y);
    end
end

%% 任意分布
Iteration2 = 32;
sum_rate_PDM_rand=zeros(Iteration2,1);
sum_rate_MF = zeros(Iteration2,1);
parfor a=1:Iteration2

    sprintf('任意位置，第%d轮',a)

    user_position_rand=zeros(3,K);       %用户坐标

    L=2;   %用户中心的Z轴坐标
    DD = 30*rand(K,1);
    angle_rand = 2*pi*rand(K,1);
    L_rand = 30*rand(K,1);

    for k = 1:K
        user_position_rand(:,k)=[DD(k)*cos(angle_rand(k)) DD(k)*sin(angle_rand(k)) L_rand(k)];
    end

    theta_k=zeros(3,K,NumdS);
    temp=0;
    for n=1:NumdS
        for k=1:K
            theta_k(:,k,n)=(2*rand(3,1)-1+2j*(rand(3,1)-0.5))/sqrt(6);
            temp=temp+norm(theta_k(:,k,n))^2*dS;
        end
    end
    theta_k=theta_k*sqrt(P_T)/sqrt(temp);

    G_k = Freespace_G_k_generate(lambda,user_position_rand,Vs,K,NumdS,delta_x,delta_y);
    [~,~,sum_rate_PDM_rand(a)] = H_MIMO_precoding(K,NumdS,NumF,dS,G_k,Psi_k,theta_k,f_base,P_T,sigma_2,Ls_x,Ls_y);
    [~,~,sum_rate_MF(a)] = MF_MIMO_precoding(K,NumdS,dS,G_k,Psi_k,theta_k,P_T,sigma_2);

end

sum_rate_PDM_rand = ones(length_L,1)*mean(sum_rate_PDM_rand);
sum_rate_MF = ones(length_L,1)*mean(sum_rate_MF);

toc
%%
sum_rate_PDM_1=mean(sum_rate_PDM_1,2);
sum_rate_PDM_2=mean(sum_rate_PDM_2,2);
sum_rate_PDM_3=mean(sum_rate_PDM_3,2);
sum_rate_PDM_4=mean(sum_rate_PDM_4,2);

sum_rate_MF_1 = mean(sum_rate_MF_1,2);
% angle = angle/pi*180;

figure;
hold on;
grid on;
box on;

Buchang = 1;

a = plot(distance(1:Buchang:end),sum_rate_PDM_1(1:Buchang:end),'-rp','LineWidth',1.5);
b = plot(distance(1:Buchang:end),sum_rate_PDM_2(1:Buchang:end),':^','LineWidth',1.5,'color',String2Color('#ffa500')/256);
c = plot(distance(1:Buchang:end),sum_rate_PDM_3(1:Buchang:end),'-.+','LineWidth',1.5,'color',String2Color('#32cd32')/256);
d = plot(distance(1:Buchang:end),sum_rate_PDM_4(1:Buchang:end),'-x','LineWidth',1.5,'color',String2Color('#8b0000')/256);
e = plot(distance,sum_rate_PDM_rand,'-r','LineWidth',1.5);

aa = plot(distance(1:Buchang:end),sum_rate_MF_1(1:Buchang:end),'-bs','LineWidth',1.5);
f = plot(distance,sum_rate_MF,':b','LineWidth',1.5);

xlabel('Network radius $R$ (m)','Interpreter','latex');
ylabel('Sum-rate (bps/Hz)','Interpreter','latex');
legend([a b c d e aa f],'PDM $(L = 2~{\rm m})$', 'PDM $(L = 5~{\rm m})$','PDM $(L = 10~{\rm m})$','PDM $(L = 30~{\rm m})$','PDM (Random positions)','MF $(L = 2~{\rm m})$','MF (Random positions)','Interpreter','latex','FontSize',12);
set(gca,'FontName','Times','FontSize',14);

