% WEBDEMO .m-file (Octave), Institute of Telecommunications, University of Stuttgart, Germany
% AUTHOR:   Nan Zhao

function [] = sinr(
    diagram,                    % ccdf = 1, cdf = 2, pdf = 3
    nChannel,                   % number of channel
    nSymbol,                    % number of symbol
    snr,                        % in dB, start value of SNR
    nTx,                        % number of transmit antenna
    nRx,                        % number of receive antenna
    ZF_en,                      % enable ZF
    MMSE_en,                    % enable MMSE
    ZF_SIC_en,                	% enable ZF with unordered SIC
    MMSE_SIC_en,              	% enable MMSE with  SIC
    snr_range,                  % ± 10dB
    snrCompare_en,              % observe the curve by comparing two neigboring snr
    svg_en,                     % if 1, then also output svg-file
    webdemo_en,                 % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
    filename                    % output file for webdemo
)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN, set up graphical output

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

    landscape = "-S930,350";
    portrait  = "-S640,480";

    output_format = portrait;     % default

    set(0, "defaultlinelinewidth", 2);

      if (webdemo_en)
        set(0, "defaultaxesfontsize", 8);
        set(0, "defaultaxesfontname", "Helvetica");
        set(0, "defaulttextfontsize", 8);
        set(0, "defaulttextfontname", "Helvetica");
        const_lw = 1;
        const_msize = 4;
      elseif (svg_en)
        set(0, "defaultaxesfontsize", 6);
        set(0, "defaultaxesfontname", "Helvetica");
        set(0, "defaulttextfontsize", 6);
        set(0, "defaulttextfontname", "Helvetica");
        const_lw = 1;
        const_msize = 3;
      else
        set(0, "defaultaxesfontsize", 10);
        set(0, "defaultaxesfontname", "Helvetica");
        set(0, "defaulttextfontsize", 8);
        set(0, "defaulttextfontname", "Helvetica");
        const_lw = 2;
        const_msize = 2;
      endif

    % Set default color order for multiline plots
    set(0, 'DefaultAxesColorOrder', unicmap);

    % END, set up graphical output
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN check some input parameters for consistency

    userWarning = 0;
    warnStr=['Warning: ',"\n"];
    if(nChannel>1000000)
        nChannel = 100000;
        %error('too much channel.');
        warnStr= [warnStr "reset the number of channel to 100000\n"];
        userWarning=1;
    endif

    if(nRx<nTx)
        nRx = nTx;
        %error('receive antennas cannot be smaller than transmit antennas');
        warnStr= [warnStr "reset the number of receive antennas equal to the number of transmit antennas\n"];
        userWarning=1;
    endif

    if(nTx>4)
        if(ZF_SIC_en || MMSE_SIC_en)
            %error('too much transmit antennas');
            ZF_SIC_en = 0;
            MMSE_SIC_en = 0;
            ZF_en = 1;
            MMSE_en = 1;
            warnStr= [warnStr "ML and SIC methode donot support over 4 transmit antennas\n"];
            userWarning=1;
        endif
    endif

    % END
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    clf
    figure(1);
    hold on;
    grid on;

    if(!snrCompare_en)
        mainFunction(snr,snr,diagram,nChannel,nTx,nRx,nSymbol,ZF_en,MMSE_en,ZF_SIC_en,MMSE_SIC_en,snr_range,1);
    endif

    if(snrCompare_en)
        mainFunction(snr,snr,diagram,nChannel,nTx,nRx,nSymbol,ZF_en,MMSE_en,ZF_SIC_en,MMSE_SIC_en,snr_range,1);
        mainFunction(snr+10,snr,diagram,nChannel,nTx,nRx,nSymbol,ZF_en,MMSE_en,ZF_SIC_en,MMSE_SIC_en,snr_range,0);
        mainFunction(snr-10,snr,diagram,nChannel,nTx,nRx,nSymbol,ZF_en,MMSE_en,ZF_SIC_en,MMSE_SIC_en,snr_range,0);
    endif

    % Warning message of invalid parameter settings
    if (userWarning)
     text(0.05, 0.5, warnStr, 'Color', 'red', 'FontSize', 7,...
            'Units', 'Normalized', 'HorizontalAlignment', 'left');
    end 
    

    % Create an image file. This image is displayed on the webpage. Its size
    % must be 640x480.
    %print(filename, "-dgif","-S640,480", "-r0");
    print(filename, "-dpng",output_format);
    %print(filename, "-djpg","-S1280,960", "-r0");
    if (svg_en)
    print(["fig_",filename], "-dsvg",output_format);
    endif

endfunction

% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main function

function [] = mainFunction(snr,snr_original,diagram,nChannel,nTx,nRx,nSymbol,ZF_en,MMSE_en,ZF_SIC_en,MMSE_SIC_en,snr_range,legend_en)
    
    % pre-allocate memory for outputs
    nSinrZF = zeros(nChannel,nTx);
    nSinrMMSE = zeros(nChannel,nTx);
    nSinrZF_UOSIC = zeros(nChannel,nTx);
    nSinrZF_OSIC = zeros(nChannel,nTx);
    nSinrMMSE_UOSIC = zeros(nChannel,nTx);
    nSinrMMSE_OSIC = zeros(nChannel,nTx);

    step = 1;
    range = -snr_range:step:snr_range;
    plotstreamIdx = 1;

    % init parameters for plot

    if(diagram == 1)
        ylabelString = 'CCDF';
        axis([-snr_range snr_range,0,1.1]);
        legendLocation = 'SouthWest';
    end

    if(diagram == 2)
        ylabelString = 'CDF';
        axis([-snr_range snr_range,0,1.1]);
        legendLocation = 'SouthEast';
    end

    if(diagram == 3)
        ylabelString = 'PDF';
        axis([-snr_range snr_range,0,0.15]);
        legendLocation = 'NorthWest';
    end

    xlabel('SINR in dB');
    ylabel(ylabelString);
    title(ylabelString);

    init_legend = ['1';'2';'3';'4';'5';'6'];
    legend_str = cellstr(init_legend);
    label_legend = cellstr(init_legend);

    % main modul
    for ii = 1:nChannel
    
        [s_tx,sigma_n,h,r] = generate_MIMO_signal(snr,nTx,nRx,nSymbol);
    
        if(ZF_en)
            nSinrZF(ii,:) = sinrZF(h,nTx,sigma_n);
        end

        if(MMSE_en)
            nSinrMMSE(ii,:) = sinrMMSE(h,nTx,sigma_n);
        end
        
        if(ZF_SIC_en)
            nSinrZF_UOSIC(ii,:) = sinrSIC(s_tx,h,r,nTx,sigma_n,0,0);
            nSinrZF_OSIC(ii,:) = sinrSIC(s_tx,h,r,nTx,sigma_n,0,1);
        end

        if(MMSE_SIC_en)
            nSinrMMSE_UOSIC(ii,:) = sinrSIC(s_tx,h,r,nTx,sigma_n,1,0);
            nSinrMMSE_OSIC(ii,:) = sinrSIC(s_tx,h,r,nTx,sigma_n,1,1);
        end
    end

    % plot
    if(ZF_en)
        outputZF = get_distribution(10*log10(nSinrZF(:,plotstreamIdx)),diagram,range);
        plot(range,outputZF,"color",unicmap(1));
        legend_str(1) = cellstr('ZF');
    else
        legend_str(1) = cellstr('');
    end

    if(MMSE_en)
        outputMMSE = get_distribution(10*log10(nSinrMMSE(:,plotstreamIdx)),diagram,range);
        plot(range,outputMMSE,"color",unicmap(2));
        legend_str(2) = cellstr('MMSE');
    else
        legend_str(2) = cellstr('');
    end

    if(ZF_SIC_en)
        outputZF_UOSIC = get_distribution(10*log10(nSinrZF_UOSIC(:,plotstreamIdx)),diagram,range);
        outputZF_OSIC = get_distribution(10*log10(nSinrZF_OSIC(:,plotstreamIdx)),diagram,range);

        plot(range,outputZF_UOSIC,"-:","color",unicmap(3));
        plot(range,outputZF_OSIC,"-:","color",unicmap(7));
        legend_str(3) = cellstr('ZF with unordered SIC');
        legend_str(4) = cellstr('ZF with ordered SIC');

    else
        legend_str(3) = cellstr('');
        legend_str(4) = cellstr('');
    end

    if(MMSE_SIC_en)
        outputMMSE_UOSIC = get_distribution(10*log10(nSinrMMSE_UOSIC(:,plotstreamIdx)),diagram,range);
        outputMMSE_OSIC = get_distribution(10*log10(nSinrMMSE_OSIC(:,plotstreamIdx)),diagram,range);

        plot(range,outputMMSE_UOSIC,"--","color",unicmap(5));
        plot(range,outputMMSE_OSIC,"--","color",unicmap(6));
        legend_str(5) = cellstr('MMSE with unordered SIC');
        legend_str(6) = cellstr('MMSE with ordered SIC');
    else
        legend_str(5) = cellstr('');
        legend_str(6) = cellstr('');
    end

    

    idx=1;
    for i = 1:6
        if(~isempty(char(legend_str(i))))
            label_legend(idx) = legend_str(i);
            idx = idx+1;
        end
    end
    
    if(legend_en)
        legend({char(label_legend(1)),char(label_legend(2)),char(label_legend(3)),char(label_legend(4)),char(label_legend(5)),char(label_legend(6))},'Location',legendLocation);      
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% channel simulator
    function [ s_tx,sigma_n,h,y ] = generate_MIMO_signal(snr,nTx,nRx,nSymbol)
        sigma_n = 10^(-snr/10);
        s_tx = 1/sqrt(2)*(sign(randn(nTx,nSymbol))+1i*(sign(randn(nTx,nSymbol))))
        h = (normrnd(0,1/sqrt(2),nRx,nTx) + 1i.*normrnd(0,1/sqrt(2),nRx,nTx))/sqrt(nTx);
        n = normrnd(0,sqrt(sigma_n./2),nRx,nSymbol) + 1i.*normrnd(0,sqrt(sigma_n./2),nRx,nSymbol);
        y = h*s_tx+n;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ZF
	function [nSinrZF] = sinrZF(h,nTx,sigma_n)
	    w_ZF = (h'*h)\h';

        P_ZF = ones(nTx,1);
        N_ZF = sigma_n*diag(real(w_ZF*w_ZF'));
	    nSinrZF = P_ZF./N_ZF;
	end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MMSE
	function [nSinrMMSE] = sinrMMSE(h,nTx,sigma_n)
	    w_MMSE = (h'*h+sigma_n*eye(nTx))\h';
	    d_MMSE = diag(diag(w_MMSE*h));
	    z_MMSE = w_MMSE*h-d_MMSE;
	    
	    P_MMSE = diag(real(d_MMSE*d_MMSE'));
        I_MMSE = diag(real(z_MMSE*z_MMSE'));
        N_MMSE = diag(sigma_n*real(w_MMSE*w_MMSE'));
	    
	    nSinrMMSE = P_MMSE./(I_MMSE+N_MMSE);
	end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% to detemine the range of soft value
    function [s_est_hard] = HardDecisionUnit(s_est_soft)
        s_est_hard = 1/sqrt(2)*(sign(real(s_est_soft))+1i*sign(imag(s_est_soft)));
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ nSinr ] = sinrSIC(s,h,r,nTx,sigma_n,useMMSE,activeOrdering)

    % h is the channel matrix
    % r is the received undemapped signal
    % methode means the called algorithm of detections, ZF or MMSE

    mmse_enable = useMMSE;
    nSinr = zeros(nTx,1);
    nSymbol =size(s,2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % first iteration
    % calculate the product of receive matrix W's first row and the signal
    % vector r, obtaining the soft-detected value of first layer 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    w_1 = (h'*h+sigma_n*eye(nTx)*mmse_enable)\h';
    %select strongst layer to detect
    if(activeOrdering)
        streamIdx_1 = get_first_SIC_stream_todetect(w_1,h,sigma_n);
    else
        streamIdx_1 = 1;
    end
    w1_1 = w_1(streamIdx_1,:);
    y_1 = w1_1*r;

    Noise_P = sigma_n*real(w1_1*(w1_1')); % noise power
    Sig_tmp = w1_1*h;
    Sig_P = norm(Sig_tmp(streamIdx_1))^2;
    Interference_P = norm(Sig_tmp)^2-Sig_P;
    nSinr(1) = Sig_P/(Interference_P+Noise_P);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % second iteration
    % h2 = [h2,h3,...] h2,h3 are the ranks of H
    % substracte the detected signal and its corresponding parts from the
    % channel matrix h and the signal vector r.
    % if the <activeOrdering> is actived, the channel matrix h should be sorted
    % according to the SINR of different layers.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    s_est_hard_1 = HardDecisionUnit(y_1,1);

    h1 = h(:,streamIdx_1);
    h_2 = h;
    h_2(:,streamIdx_1) = [];
    w_2 = (h_2'*h_2+sigma_n*eye(nTx-1)*mmse_enable)\h_2';
    %select strongst layer to detect
    if(activeOrdering)
        streamIdx_2 = get_first_SIC_stream_todetect(w_2,h_2,sigma_n);
    else
        streamIdx_2 = 1;
    end
    w1_2 = w_2(streamIdx_2,:);
    r_2 = r-h1*s_est_hard_1;
    y_2 = w1_2*r_2;

    Noise_P = sigma_n*real(w1_2*(w1_2')); % noise power
    Sig_tmp = w1_2*h_2;
    Sig_P = norm(Sig_tmp(streamIdx_2))^2;
    Interference_P = norm(Sig_tmp)^2-Sig_P;
    s1 = s(streamIdx_1,:);
    s(streamIdx_1,:)=[];
    Error_prop_P = norm(w1_2*h1*(s_est_hard_1-s1))^2/nSymbol;
    nSinr(2) = Sig_P/(Interference_P+Noise_P+Error_prop_P);


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % third iteration,repeat the process of 2nd iteration
    % H3 = [h3,h4...]
    % w_3 = inv(H3+sigma_n*I(2))
    % r_3 = r_2-h2*s^_2  h2 is the 2nd rank of channel matrix
    % y_3 = w1_3*r_3     w1_3 is the first row of w_3
    % s^_3 is the result of hard decision from y_4
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % nTx >= 3
    if nTx == 3 || nTx == 4
        s_est_hard_2 = HardDecisionUnit(y_2,1);
        h2 = h_2(:,streamIdx_2);
        h_3 = h_2;
        h_3(:,streamIdx_2) = [];
        w_3 = (h_3'*h_3+sigma_n*eye(nTx-2)*mmse_enable)\h_3';
        %select strongst layer to detect
        if(activeOrdering)
            streamIdx_3 = get_first_SIC_stream_todetect(w_3,h_3,sigma_n);
        else
            streamIdx_3 = 1;
        end
        w1_3 = w_3(streamIdx_3,:);
        r_3 = r_2-h2*s_est_hard_2;
        y_3 = w1_3*r_3;
        
        Noise_P = sigma_n*real(w1_3*(w1_3')); % noise power
        Sig_tmp = w1_3*h_3;
        Sig_P = norm(Sig_tmp(streamIdx_3))^2;
        Interference_P = norm(Sig_tmp)^2-Sig_P;
        s2 = s(streamIdx_2,:);
        s(streamIdx_2,:)=[];
        Error_prop_P = (norm(w1_3*h2*(s_est_hard_2-s2))^2 + norm(w1_3*h1*(s_est_hard_1-s1))^2)/nSymbol;
        nSinr(3) = Sig_P/(Interference_P+Noise_P+Error_prop_P);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % fourth iteration, simply repeat previous step
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nTx == 4
        s_est_hard_3 = HardDecisionUnit(y_3,1);
        h3 = h_3(:,streamIdx_3);
        h_4 = h_3;
        h_4(:,streamIdx_3) = [];
        w_4 = (h_4'*h_4+sigma_n*eye(nTx-3)*mmse_enable)\h_4';
        %select strongst layer to detect
        if(activeOrdering)
            streamIdx_4 = get_first_SIC_stream_todetect(w_4,h_4,sigma_n);
        else
            streamIdx_4 = 1;
        end
        w1_4 = w_4(streamIdx_4,:);
        r_4 = r_3-h3*s_est_hard_3;
        y_4 = w1_4*r_4;
        
        %evm(4) = (norm(y_4-s(4)))^2;
        Noise_P = sigma_n*real(w1_4*(w1_4')); % noise power
        Sig_tmp = w1_4*h_4;
        Sig_P = norm(Sig_tmp(streamIdx_4))^2;
        Interference_P = norm(Sig_tmp)^2-Sig_P;
        s3 = s(streamIdx_3,:);
        s(streamIdx_3,:)=[];
        Error_prop_P = (norm(w1_4*h3*(s_est_hard_3-s3))^2 + norm(w1_4*h2*(s_est_hard_2-s2))^2 + norm(w1_4*h1*(s_est_hard_1-s1))^2)/nSymbol;
        nSinr(4) = Sig_P/(Interference_P+Noise_P+Error_prop_P);
    end

    nSinr = (nSinr);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunktion of SIC algorithm
function streamIdx = get_first_SIC_stream_todetect(w_rx,channelH,noise_var)
    G = w_rx*channelH;
    P_Sig = abs(diag(G)).^2;
    P_ICI = (abs(diag((G-diag(diag(G)))*((G-diag(diag(G)))')))).^2;
    P_noise = noise_var*diag((w_rx*(w_rx')));
    SINR = P_Sig./(P_ICI+P_noise);
    [~,streamIdxSorted] = sort(SINR,'descend');
    streamIdx = streamIdxSorted(1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate the pdf, cdf and ccdf of after detection

    function [outputValue_diagram,range] = get_distribution(sinr,diagram,range)
        step = 1;
        histogram = hist(sinr,range);
        if(diagram == 3)
            outputValue_diagram = histogram/(step*sum(histogram)); % PDF
        end

        if(diagram == 2)
            outputValue_diagram = cumsum(histogram)/(sum(histogram)); % CDF
        end

        if(diagram == 1)
            outputValue_diagram = 1-cumsum(histogram)/(sum(histogram)); % CCDF
        end

    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function [ uni_colors ] = unicmap(n)
        %UNICMAP Get a matrix of colors defined by the University of Stuttgart
        %   unicmap returns a matrix of 9 colors defined by the University of 
        %           Stuttgart.
        %   unicmap(n) returns the University color with index n. n must be in
        %              [1..9].

        % uni_hellgrau = [224/255 224/255 224/255];
        % uni_mittelgrau = [128/255 128/255 128/255];
        % uni_dunkelgrau = [80/255 80/255 80/255];    

        narginchk(0, 1);

        uni_colors = [0 126/255 198/255]; % mittelblau
        uni_colors = [uni_colors; 238/255 28/255 35/255]; % rot
        uni_colors = [uni_colors; 140/255 198/255 62/255]; % apfelgruen
        uni_colors = [uni_colors; 1 221/255 0]; % gelb
        uni_colors = [uni_colors; 244/255 111/255 33/255]; % orange
        uni_colors = [uni_colors; 237/255 0 140/255]; % pink
        uni_colors = [uni_colors; 19/255 31/255 31/255]; % anthrazit
        uni_colors = [uni_colors; 0/255 170/255 173/255]; % turquoise
        uni_colors = [uni_colors; 128/255 40/255 145/255]; % lila

        if nargin == 1
            if (~isscalar(n) || ~isnumeric(n))
                error('First argument must be a numeric scalar.')
            end
            if (n < 1 || n > size(uni_colors, 1))
                error(['First argument mus be postive and less than or equal to '...
                    num2str(size(uni_colors, 1))]);
            end
            uni_colors = uni_colors(n,:);
        end

    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
