%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .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_3_function(p_1, p_2, p_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
hold on;
% set default Fontsize and name
set(gca, 'defaultlinelinewidth', 2);
set(gca, 'fontname','Helvetica');
set(gca, 'fontsize', 15);

p_1 = 10.^(p_1/10);
p_2 = 10.^(p_2/10); % dBm -> milliwatt
p_3 = 10.^(p_3/10);
r_1 = log2(1+p_1);
r_2 = log2(1+p_2); % boundry
r_3 = log2(1+p_3);
step=50;

Point_Group=[];
R_1 = 0:r_1/step:r_1;
R_2 = 0:r_2/step:r_2;
R_3 = r_3:-r_3/step:0;
for i=1:length(R_1)
    for j=1:length(R_2)
        label=0;
        for k=1:length(R_3)
            if(((R_1(i)+R_2(j)+R_3(k))<=log2(1+p_1+p_2+p_3))&&((R_1(i)+R_2(j))<=log2(1+p_1+p_2))&&((R_1(i)+R_3(k))<=log2(1+p_1+p_3))&&((R_2(j)+R_3(k))<=log2(1+p_2+p_3)))
                Point_Group=cat(1,Point_Group,[R_1(i) R_2(j) R_3(k)]);
                label=1;
                break;
            end
        end
        if label==0
           Point_Group=cat(1,Point_Group,[R_1(i) R_2(j) nan]);
        end
    end
end
% Prepare plot
mesh(R_1,R_2,reshape(Point_Group(:,3),[(step+1),(step+1)]))

R_1 = 0:r_1/step:r_1;
R_2 = r_2:-r_2/step:0;
R_3 = 0:r_3/step:r_3;
Point_Group=[];
for k=1:length(R_3)
    for i=1:length(R_1)
        label=0;
        for j=1:length(R_2)
            if(((R_1(i)+R_2(j)+R_3(k))<=log2(1+p_1+p_2+p_3))&&((R_1(i)+R_2(j))<=log2(1+p_1+p_2))&&((R_1(i)+R_3(k))<=log2(1+p_1+p_3))&&((R_2(j)+R_3(k))<=log2(1+p_2+p_3)))
                Point_Group=cat(1,Point_Group,[R_3(k) R_1(i) R_2(j) ]);
                label=1; 
                break;
            end
        end
        if label==0
           Point_Group=cat(1,Point_Group,[R_3(k) R_1(i) nan]);
        end
    end
end
% Prepare plot
[RR_3,RR_1] = meshgrid(R_3,R_1);
mesh(RR_1,reshape(Point_Group(:,3),[(step+1),(step+1)]),RR_3)

R_1 = r_1:-r_1/step:0;
R_2 = 0:r_2/step:r_2;
R_3 = 0:r_3/step:r_3;
Point_Group=[];
for k=1:length(R_3)
    for j=1:length(R_2)
        label=0;
        for i=1:length(R_1)
            if(((R_1(i)+R_2(j)+R_3(k))<=log2(1+p_1+p_2+p_3))&&((R_1(i)+R_2(j))<=log2(1+p_1+p_2))&&((R_1(i)+R_3(k))<=log2(1+p_1+p_3))&&((R_2(j)+R_3(k))<=log2(1+p_2+p_3)))
                Point_Group=cat(1,Point_Group,[ R_3(k) R_2(j) R_1(i) ]);
                label=1; 
                break;
            end
        end
        if label==0
           Point_Group=cat(1,Point_Group,[ R_3(k) R_2(j) nan]);
        end
    end
end
% Prepare plot
[RR_3,RR_2] = meshgrid(R_3,R_2);
mesh(reshape(Point_Group(:,3),[(step+1),(step+1)]),RR_2,RR_3)

grid on;
box off;
xlabel('R_1 in bits/s/Hz');
ylabel('R_2 in bits/s/Hz');
zlabel('R_3 in bits/s/Hz');
set(gca,'xtick',[0:1:r_1],'ytick',[0:1:r_2],'ztick',[0:1:r_3]);
view(view_horizon,view_vertical); % set image view % 120,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