function [wc_opt,ws_opt,Wc_opt,Ws_opt,SNR_BS_1] = Optimization_wc_ws(EdeltaS,r_max,Delta_theta,Delta_phi,lamda,theta_S,phi_S,HBR,h_RK1,h_RK2,h_RK3,h_BK1,h_BK2,h_BK3,sigma_BS,sigma_UE,Ptx,M,Nx,Ny,ome_in,wrx_in,d,accu,gamma)

Mr = 0;

Dff_S = ((Delta_phi*(1/accu))*(Delta_theta*(1/accu)));

for ind_phi = 1:accu

phi_tem = phi_S - Delta_phi/2 + Delta_phi*(ind_phi/accu);

for ind_theta = 1:accu

theta_tem = theta_S - Delta_theta/2 + Delta_theta*(ind_theta/accu);

% Array response

for i = 1:Nx

a_vec_tem1(i,1) = exp(1j*(2*pi*d/lamda)*(i-1)* (sin(theta_tem)*cos(phi_tem)) );
end
for i = 1:Ny
a_vec_tem2(i,1) = exp(1j*(2

*pi*d/lamda)

*(i-1)*(sin(theta_tem)

*sin(phi_tem)) );*

end

a_RIS_tem = kron(a_vec_tem1,a_vec_tem2);

% Matrix Mr:

mat_C = (EdeltaS(lamda^2)/(64*(pi^3)

end

a_RIS_tem = kron(a_vec_tem1,a_vec_tem2);

% Matrix Mr:

mat_C = (EdeltaS

*(r_max^2)*(sigma_BS)))…

*( HBR’*a_RIS_tem’)*diag(ome_in’)

*diag(ome_in’)*(a_RIS_tem*HBR*(wrx_in…

*wrx_in’)

*HBR’*a_RIS_tem’)*diag(ome_in)*HBR * sin(theta_tem) );

*diag(ome_in)*(a_RIS_tem```
Mr = Mr + mat_C * Dff_S;
end
```

end

% Matrix Xi

Ome_in = diag(ome_in);

Xi_1 = (HBR’*Ome_in’*h_RK1+h_BK1)*(h_RK1’*Ome_in*HBR+h_BK1’);

Xi_2 = (HBR’*Ome_in’*h_RK2+h_BK2)*(h_RK2’*Ome_in*HBR+h_BK2’);

Xi_3 = (HBR’*Ome_in’*h_RK3+h_BK3)*(h_RK3’*Ome_in*HBR+h_BK3’);

% Optimize wc ws

Mr = round(Mr,12);

Xi_1 = round(Xi_1,15);

Xi_2 = round(Xi_2,16);

Xi_3 = round(Xi_3,16);

cvx_begin quiet

cvx_solver sedumi

cvx_precision best

variable Wc(M,M) hermitian;

variable Ws(M,M) hermitian;

% variable wc(M,1);

% variable ws(M,1);

maximize (trace((Mr)*(Wc+Ws)) )
subject to
Wc == hermitian_semidefinite(M);
Ws == hermitian_semidefinite(M);
trace(Wc + Ws) <= Ptx;
trace((Xi_1)*(Wc - gamma

*Ws)) >= gamma*(sigma_UE);

trace((Xi_2)

*(Wc - gamma*Ws)) >= gamma*(sigma_UE);

trace((Xi_3)

*(Wc - gamma*Ws)) >= gamma*(sigma_UE);

cvx_end

Wc_opt = Wc;

Ws_opt = Ws;

% Eigenvalue Decomposition

[Wc_x,Wc_y] = eig(Wc);

lamda_vec_Wc = diag(Wc_y);

[lamda_max_Wc,ii] = max(lamda_vec_Wc);

eigvec_max_c = Wc_x(:,ii);

wc_opt = sqrt(lamda_max_Wc)*eigvec_max_c;

[Ws_x,Ws_y] = eig(Ws);

lamda_vec_Ws = diag(Ws_y);

[lamda_max_Ws,ii] = max(lamda_vec_Ws);

eigvec_max_s = Ws_x(:,ii);

ws_opt = sqrt(lamda_max_Ws)*eigvec_max_s;

SNR_BS_1 = abs(trace((Mr)*(Wc_opt + Ws_opt)) );

end