%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo 'Multiuser Capacity' 
% Written for Matlab
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Wei Wang, Jintong JU   
% date: 14.12.2020 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example script just prints the values of the given parameters.
% Make sure that the function prototype matches the prototype defined
% in the XML config file. (Parameter names are arbitrary but should be the same
% for better readability)
function [] = demo_2_function_2(p_1,p_2,sigma_min,sigma_max,view_horizon,view_vertical, webdemo_en,filename)

% Create an invisible figure.
if (webdemo_en)
    fig = figure(1); set(fig, 'visible', 'off');
else
    fig = figure(1); set(fig, 'visible', 'on');
end

% set default Fontsize and name
set(gca, 'defaultlinelinewidth', 2);
set(gca, 'fontname','Helvetica');
set(gca, 'fontsize', 20);

sig_min = 10.^(sigma_min/10);
sig_max = 10.^(sigma_max/10);
P_1 = 10.^(p_1/10);
P_2 = 10.^(p_2/10); % dBm -> milliwatt
step=(sig_min-sig_max)/10;
sigma=sig_max:step:sig_min;

R_1_array=[];
R_2_array=[];
SNR_array=[];
for i=1:length(sigma)
    r_1_temp = log2(1+P_1/sigma(i));
    r_2_temp = log2(1+P_2/sigma(i));
    r_1_left_temp = log2(1+(P_1+P_2)/sigma(i))-r_2_temp;
    R_1 = 0:r_1_temp/100:r_1_temp;
    R_2 = r_2_temp.*(R_1<r_1_left_temp)+(log2(1+(P_1+P_2)/sigma(i))-R_1).*(R_1>=r_1_left_temp&R_1<r_1_temp)+0.*(R_1>=r_1_temp);
    R_1_array=cat(1,R_1_array,R_1);
    R_2_array=cat(1,R_2_array,R_2);
    snr_temp=10*log10(P_1/sigma(i));
    SNR_array=cat(1,SNR_array,snr_temp*ones(1,101)); 
end

surf(SNR_array,R_1_array,R_2_array);
xlabel({'SNR of user 1';'       in dB'},'fontsize',15);
ylabel('R_1 in bits/s/Hz','fontsize',15);
zlabel('R_2 in bits/s/Hz','fontsize',15);

set(gca,'YDir','reverse');
colorbar('fontname','Helvetica','fontsize',12);
title(colorbar,'rate [bits/s/Hz]','fontname','Helvetica','fontsize',12);
set(gca,'xtick',[SNR_array(1):2:SNR_array(11)],'ytick',[0:1:r_1_temp],'ztick',[0:1:r_2_temp]);
view(view_horizon,view_vertical); % set image view % -60,20
axis equal;

% axesLabelsAlign3D;
% set the lables parallel to the axes
% copyright Matthew Arthington
[az,el] = view;
Raz = [cosd(az) sind(az) 0;...
	-sind(az) cosd(az) 0;...
	0 0 1];
Rel = [1 0 0;...
	0 cosd(el) -sind(el);
	0 sind(el) cosd(el)];
u=get(gca,'CameraUpVector');
if ~(u(1)==0 && u(2)==0 && u(3)==1)
	p=get(gca,'CameraPosition');
	t=get(gca,'CameraTarget');
    v=(p-t);
    v=v/norm(v);%View vector from camera to target
    u = u/norm(u);%Camera up vector
    q = cross(v,u);
    q=q/norm(q);
    %Get the x axis's projection into the view plane and then find its angle wrt the up vector
    xH = cross([1;0;0],v);
    xH = xH/norm(xH);
    xH = cross(xH,v);
    thetax = -acosd(dot(xH,u))+90;
    %Check which way the label needs to be rotated
    if dot(q,xH)>0
        thetax = -thetax;
    end
    %Get the y axis's projection into the view plane and then find its angle wrt the up vector
    yH = cross([0;1;0],v);
    yH = yH/norm(yH);
    yH = cross(yH,v);
    thetay = -acosd(dot(yH,u))+90;
    if dot(q,yH)>0
        thetay = -thetay;
    end
    %Get the z axis's projection into the view plane and then find its angle wrt the up vector
    zH = cross([0;0;1],v);
    zH = zH/norm(zH);
    zH = cross(zH,v);
    thetaz = -acosd(dot(zH,u))+90;
    if dot(q,zH)>0
        thetaz = -thetaz;
    end
else %When rotate3d has been used, the up vector isn't set by matlab correctly
    %Calculate current orientation of x and y axes in view coordinates
    xax = Rel*Raz*[1;0;0];
    yax = Rel*Raz*[0;1;0];
    %Project x and y into current viewing plane
    n1=cross(xax,[0;1;0]);
    x = cross([0;1;0],n1);
    n1=cross(yax,[0;1;0]);
    y = cross([0;1;0],n1);
    thetax = atand(x(3)/x(1));
    thetay = atand(y(3)/y(1));
    if ~any(x)
        thetax = 0;
    end
    if ~any(y)
        thetay = 0;
    end
    thetaz = 90;
end
%Orientate these labels to be aligned with the
%axis directions.
set(get(gca,'xlabel'),'rotation',thetax);
set(get(gca,'ylabel'),'rotation',thetay);
set(get(gca,'zlabel'),'rotation',thetaz);

axis normal;

if (webdemo_en)
    print(filename, "-dpng",'-r85');% generate the picture results
end

end