clear all; clc;

addpath(genpath('H-MIMO-functions'));
addpath(genpath('Upperbound-functions'));

tic

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

Z0=376.73;              %真空阻抗
sigma_2=5.6*10^(-3)/kappa^2/Z0^2;    %V2/m2
P_T=10^(-4);                        %电流密度功率P_T

K=8;                    %用户数量

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

L=30;   %用户中心的Z轴坐标

user_position(:,1)=[1 1 L];
user_position(:,2)=[1 -1 L];
user_position(:,3)=[-1 -1 L];
user_position(:,4)=[-1 1 L]; 
user_position(:,5)=[5 5 L];
user_position(:,6)=[5 -5 L];
user_position(:,7)=[-5 -5 L];
user_position(:,8)=[-5 5 L]; 

%% 发射机面积 坐标
% Ls_m=[0.2:0.2:1.6];

% L_m=[0:0.1:0.9];
% Ls_Mianji=L_m.^2;

Ls_Mianji=[0:0.1:1];
Ls_Mianji(1)=0.001;

length_L=length(Ls_Mianji);
delta_x=33*ones(length_L,1);

Iteration=1;

sum_rate_H_MIMO_Upper=zeros(length_L,Iteration);
sum_rate_H_MIMO_freespace=zeros(length_L,Iteration);
sum_rate_H_MIMO_wavenumber=zeros(length_L,Iteration);

Ls_x=zeros(length_L,1);
Ls_y=zeros(length_L,1);

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

Ls_x(m)=sqrt(Ls_Mianji(m));                           %发射机板面宽度
Ls_y(m)=sqrt(Ls_Mianji(m));
%Ls_y(m)=0.2;

delta_y=delta_x;                         %y轴量化
s_center= [0 0 0]';                 %发射机中心

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

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

%% 傅里叶基展开
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(m)*Ls_y(m))*exp(1j*2*pi*(u-(delta_u/2))*(Vs(1,n)-Ls_x(m)/2)/Ls_x(m))*exp(1j*2*pi*(v-(delta_v/2))*(Vs(2,n)-Ls_y(m)/2)/Ls_y(m));      
        end  
    end
end

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);

%% 接收机面积 坐标
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

%% 容量计算
G_k = Freespace_G_k_generate(lambda,user_position,Vs,K,NumdS,delta_x(m),delta_y(m));

[~,~,sum_rate_H_MIMO_Upper(m,q)] = Upperbound_H_MIMO_precoding2(K,NumdS,NumF,dS,G_k,Psi_k_opt,theta_k_opt,f_base,P_T,sigma_2,Ls_x(m),Ls_y(m))

end
end
toc
%%
sum_rate_H_MIMO_freespace=mean(sum_rate_H_MIMO_freespace,2);
sum_rate_H_MIMO_wavenumber=mean(sum_rate_H_MIMO_wavenumber,2);
sum_rate_H_MIMO_Upper=mean(sum_rate_H_MIMO_Upper,2);

%Ls_Mianji=Ls_Mianji*10000;

save('main_3.mat','Ls_Mianji','Ls_x','sum_rate_H_MIMO_Upper','sum_rate_H_MIMO_freespace','sum_rate_H_MIMO_wavenumber');

figure;
hold on;
grid on;
box on;
plot(Ls_x,sum_rate_H_MIMO_Upper,'--k','LineWidth',1.5);
plot(Ls_x,sum_rate_H_MIMO_freespace,'-rp','LineWidth',1.5);
plot(Ls_x,sum_rate_H_MIMO_wavenumber,'-bs','LineWidth',1.5);
xlabel('Aperture length $L_x$ and width $L_y$ ($\rm m$)','Interpreter','latex');
ylabel('Capacity (bps/Hz)','Interpreter','latex');
legend('Upper bound (without interference)','Proposed EM-waveforming','Wavenumber-division multiplexing','Interpreter','latex','FontSize',12);
set(gca,'FontName','Times','FontSize',14);
%xlim([0.1 1])
%ylim([0 25]);
%Ls_Mianji=Ls_Mianji/10000;

Buchang=1;
figure;
hold on;
grid on;
box on;
plot(Ls_Mianji,sum_rate_H_MIMO_Upper,'--k','LineWidth',1.5);
plot(Ls_Mianji,sum_rate_H_MIMO_freespace,'-r','LineWidth',1.5);
plot(Ls_Mianji,sum_rate_H_MIMO_wavenumber,'-b','LineWidth',1.5);
plot(Ls_Mianji(1:Buchang:end),sum_rate_H_MIMO_freespace(1:Buchang:end),'pr','LineWidth',1.5);
plot(Ls_Mianji(1:Buchang:end),sum_rate_H_MIMO_wavenumber(1:Buchang:end),'bs','LineWidth',1.5);

a=plot(Ls_Mianji(end),sum_rate_H_MIMO_Upper(end),'--k','LineWidth',1.5);
b=plot(Ls_Mianji(end),sum_rate_H_MIMO_freespace(end),'-pr','LineWidth',1.5);
c=plot(Ls_Mianji(end),sum_rate_H_MIMO_wavenumber(end),'-bs','LineWidth',1.5);

% set(gca,'XTick',[0:0.1:1]);
% xlim([0 1])
% ylim([0 5])
xlabel('Aperture area $A_{\rm T}$ (${\rm m}^2$)','Interpreter','latex');
ylabel('Capacity (bps/Hz)','Interpreter','latex');
legend([a b c],'Upper bound (without interference)','Pattern-division multiplexing (proposed)','Wavenumber-division multiplexing [14]','Interpreter','latex','FontSize',12);
set(gca,'FontName','Times','FontSize',14);
