%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .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_6_function(p,h_1,h_2,h_3,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', 14);

p = 10.^(p/10);
N_0=1; % milliwatt

R_1_NOMA=[];
R_2_NOMA=[];
R_3_NOMA=[];
step_3=0.01;
step_2_num=21;
beta_3 = 0:step_3:1;
for i=1:(length(beta_3)-1)
    r_3_noma=log2(1+p*beta_3(i)*h_3/N_0);
    step_2=(1-beta_3(i))/(step_2_num-1);
    beta_2=0:step_2:(1-beta_3(i));
    beta_1=(1-beta_3(i))*ones(1,step_2_num)-beta_2;
    r_1_noma=log2(1+(p.*beta_1*h_1)./(p.*(beta_2+beta_3(i))*h_1+N_0));
    r_2_noma=log2(1+(p.*beta_2*h_2)./(p*beta_3(i)*h_2+N_0));
    R_1_NOMA=cat(1,R_1_NOMA,r_1_noma);
    R_2_NOMA=cat(1,R_2_NOMA,r_2_noma);
    R_3_NOMA=cat(1,R_3_NOMA,r_3_noma.*ones(1,step_2_num));
end
% Plot into the figure
surf(R_1_NOMA,R_2_NOMA,R_3_NOMA);
view(view_horizon,view_vertical);

% set plot
xlabel('R_1 in bits/s/Hz','fontsize', 15);
ylabel('R_2 in bits/s/Hz','fontsize', 15);
zlabel('R_3 in bits/s/Hz','fontsize', 15);
set(gca,'YDir','reverse'); % set direction of Y axis
grid on;
hold on;
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