%Student: Santiago Pagola
%LiU-ID & Personal Number: 930302-T092, sanpa993
 
%IMPORTANT NOTE: AN EXTERNAL FUNCTION WAS CREATED IN ORDER TO WRAP THE SUMS
%OF THE VECTORS: a = sumk(vector,k). It will be described in the report,
%but basically sums the first k elements of an array called 'vector'.
 
 
%clear all variables from workspace, close all figures and clear everything
%in command line
clear all
close all
clc
 
%--------------------------CONSTANT DEFINITIONS-----------------------------
m = 100; %100 nodes on no-buffering assumption
t = 1000; %Time: 1000 slots of time
slots = 1:t; %Slots array
backlog_array = zeros(size(1:t)); %Array of the values of the system backlog
n = 0; %Initially, all nodes are unbacklogged
packets_leaving = 1:t; %Array which will count #packets leaving the system
packets_arriving = 1:t; %Array which will count #packets entering the system
success = 0; %number of succesful transmissions
ps_simulation = 0; %Simulated probability of success
ps_theory = 0; %Theoretical probability of success
state_probs = zeros(size(1:m)); %Array with the values of the steady state probabilities
states = 1:m; %States array
Att_Rate = 1:m; %Attemp-rate array
G_n = zeros(size(1:m)); %Attempt rate array
Ps = zeros(size(1:m));
%---------------------------------------------------------------------------
 
qr = input('Type retransmission probability, qr: '); %Ask user for qr
lambda = input('Type overal arrival rate: '); %Ask user for overall arrival rate
qa = 1-exp(1)^(-lambda/m); %qa is the probability that each unbacklogged node has to transmit in the next slot
 
 
for i = 1:t
    %Let's now calculate the probabilities that up to 10 unbacklogged nodes
    %transmit (up to 10 new arrivals)(Same with the probabilities that up to 10 backlogged nodes
    %retransmit (up to 10 new arrivals))
    Qa = zeros(size(1,11));
    Qr = zeros(size(1,11));
    for j = 0:10
        Qa (j+1) = binopdf(j,m-n,qa);
    end
    for j = 0:10
        Qr (j+1) = binopdf(j,n,qr);
    end
 
    %a and b will act as two random numbers between 0 and 1, which will be
    %two different realizations of Qa_i and Qr_i
    a = rand(1);
    b = rand(1);
 
    %Now, two conditions must be made: one for the case when there are
    %backlogged nodes, and other one for the case when there are no
    %backlogged nodes:
 
    %CASE 1: NO BACKLOGGED NODES
    if n == 0
        if 0 <= a && a <= sumk(Qa,1) %CASE: NO BACKLOGGED NODES NOR UNBACKLOGGED PACKETS, NO TRANSMISSION->IDLE SLOT
            packets_arriving(i)=0;
            packets_leaving(i)=0;
        elseif sumk(Qa,1) < a && a <= sumk(Qa,2) %CASE: NO BACKLOGGED NODES % 1 UNBACKLOGGED PACKET TRANSMITS SUCCESFULLY
            packets_arriving(i)=1;
            packets_leaving(i)=1;
            success=success+1;
        elseif sumk(Qa,2) < a && a <= sumk(Qa,3) %CASE: NO BACKLOGGED NODES & 2 UNBACKLOGGED PACKETS TRANSMIT->COLLISION
            n=n+2;
            packets_arriving(i)=2;
            packets_leaving(i)=0;
        elseif sumk(Qa,3) < a && a <= sumk(Qa,4)  %COLLISION
            n=n+3;
            packets_arriving(i)=3;
            packets_leaving(i)=0;
        elseif sumk(Qa,4) < a && a <= sumk(Qa,5) %COLLISION
            n=n+4;
            packets_arriving(i)=4;
            packets_leaving(i)=0;
        elseif sumk(Qa,5) < a && a <= sumk(Qa,6) %AND SO ON...
            n=n+5;
            packets_arriving(i)=5;
            packets_leaving(i)=0;
        elseif sumk(Qa,6) < a && a <= sumk(Qa,7)
            n=n+6;
            packets_arriving(i)=6;
            packets_leaving(i)=0;
        elseif sumk(Qa,7) < a && a <= sumk(Qa,8)
            n=n+7;
            packets_arriving(i)=7;
            packets_leaving(i)=0;
        elseif sumk(Qa,8) < a && a <= sumk(Qa,9)
            n=n+8;
            packets_arriving(i)=8;
            packets_leaving(i)=0;
        elseif sumk(Qa,9) < a && a <= sumk(Qa,10)
            n=n+9;
            packets_arriving(i)=9;
            packets_leaving(i)=0;
        elseif sumk(Qa,10) < a && a <= 1
            n=n+10;
            packets_arriving(i)=10;
            packets_leaving(i)=0;
        end   
 
    else %CASE 2: BACKLOGGED NODES. INSIDE THIS CASE, THREE OPTIONS:
        %2.1.NO RETRANSMISSION AT ALL
        %2.2.RETRANSMISSION OF 1 PACKET
        %2.3.RETRANSMISSION OF 2 OR MORE PACKETS
 
        if 0 <= a && a <= sumk(Qr,1)    % CASE 2.1: NO RETRANSMISSION OF ANY BACKLOGGED PACKETS
 
            if 0 <= b && b <= sumk(Qa,1)   %IDLE SLOT, NO NEW ARRIVALS
                packets_arriving(i)=0;
                packets_leaving(i)=0;
            elseif sumk(Qa,1) < b && b <= sumk(Qa,2) %1 NEW PACKET SUCCESFULLY ARRIVED
                packets_arriving(i)=1;
                packets_leaving(i)=1;
                success=success+1;
            elseif sumk(Qa,2) < b && b <= sumk(Qa,3) %2 NEW PACKETS ARRIVED - COLLISION
                n=n+2;
                packets_arriving(i)=2;
                packets_leaving(i)=0;
            elseif sumk(Qa,3) < b && b <= sumk(Qa,4) %3 NEW PACKETS ARRIVED - COLLISION
                n=n+3;
                packets_arriving(i)=3;
                packets_leaving(i)=0;
            elseif sumk(Qa,4) < b && b <= sumk(Qa,5) %AND SO ON...
                n=n+4;
                packets_arriving(i)=4;
                packets_leaving(i)=0;
            elseif sumk(Qa,5) < b && b <= sumk(Qa,6)
                n=n+5;
                packets_arriving(i)=5;
                packets_leaving(i)=0;
            elseif sumk(Qa,6) < b && b <= sumk(Qa,7)
                n=n+6;
                packets_arriving(i)=6;
                packets_leaving(i)=0;
            elseif sumk(Qa,7) < b && b <= sumk(Qa,8)
                n=n+7;
                packets_arriving(i)=7;
                packets_leaving(i)=0;
            elseif sumk(Qa,8) < b && b <= sumk(Qa,9)
                n=n+8;
                packets_arriving(i)=8;
                packets_leaving(i)=0;
            elseif sumk(Qa,9) < b && b <= sumk(Qa,10)
                n=n+9;
                packets_arriving(i)=9;
                packets_leaving(i)=0;
            elseif sumk(Qa,10) < b && b <= 1
                n=n+10;
                packets_arriving(i)=10;
                packets_leaving(i)=0;
            end            
 
 
        elseif sumk(Qr,1) < a && a <= sumk(Qr,2)    % CASE 2.2: RETRANSMISSION OF A SINGLE PACKET
 
            if 0 <= b && b <= sumk(Qa,1) %CASE: SUCCESFUL BACKLOG
                n=n-1; %n DECREASES 1 UNIT
                packets_arriving(i)=1;
                packets_leaving(i)=1;
                success=success+1;
            elseif sumk(Qa,1) < b && b <= sumk(Qa,2) %CASE: COLLISION->1 NEW ARRIVAL, 1 BACKLOG
                n=n+1;
                packets_arriving(i)=2;
                packets_leaving(i)=0;
            elseif sumk(Qa,2) < b && b <= sumk(Qa,3) %CASE: COLLISION->2 NEW ARRIVALS, 1 BACKLOG
                n=n+2;
                packets_arriving(i)=3;
                packets_leaving(i)=0;
            elseif sumk(Qa,3) < b && b <= sumk(Qa,4) %AND SO ON...
                n=n+3;
                packets_arriving(i)=4;
                packets_leaving(i)=0;
            elseif sumk(Qa,4) < b && b <= sumk(Qa,5)
                n=n+4;
                packets_arriving(i)=5;
                packets_leaving(i)=0;
            elseif sumk(Qa,5) < b && b <= sumk(Qa,6)
                n=n+5;
                packets_arriving(i)=6;
                packets_leaving(i)=0;
            elseif sumk(Qa,6) < b && b <= sumk(Qa,7)
                n=n+6;
                packets_arriving(i)=7;
                packets_leaving(i)=0;
            elseif sumk(Qa,7) < b && b <= sumk(Qa,8)
                n=n+7;
                packets_arriving(i)=8;
                packets_leaving(i)=0;
            elseif sumk(Qa,8) < b && b <= sumk(Qa,9)
                n=n+8;
                packets_arriving(i)=9;
                packets_leaving(i)=0;
            elseif sumk(Qa,9) < b && b <= sumk(Qa,10)
                n=n+9;
                packets_arriving(i)=10;
                packets_leaving(i)=0;
            elseif sumk(Qa,10) < b && b <= sumk(Qa,11)
                n=n+10;
                packets_arriving(i)=11;
                packets_leaving(i)=0;
            end 
 
 
        elseif sumk(Qr,2) < a && a <= 1 %CASE 2.3: RETRANSMISSION OF 2 OR MORE PACKETS
 
            if sumk(Qr,2) < a && a <= sumk(Qr,3) %Let aux be an auxiliary variable indicating the number of backlogged nodes.
                aux=2;
            elseif sumk(Qr,3) < a && a <= sumk(Qr,4)
                aux=3;
            elseif sumk(Qr,4) < a && a <= sumk(Qr,5)
                aux=4;
            elseif sumk(Qr,5) < a && a <= sumk(Qr,6)
                aux=5;
            elseif sumk(Qr,6) < a && a <= sumk(Qr,7)
                aux=6;
            elseif sumk(Qr,7) < a && a <= sumk(Qr,8)
                aux=7;
            elseif sumk(Qr,8) < a && a <= sumk(Qr,9)
                aux=8;
            elseif sumk(Qr,9) < a && a <= sumk(Qr,10)
                aux=9;
            end
 
            if  0 <= b && b <= sumk(Qa,1)  % 0 ARRIVALS->BACKLOG REMAINS THE SAME BEACUSE NO UNBACKLOGGED NODES TRANSMIT 
                packets_arriving(i)=aux;
                packets_leaving(i)=0;
            elseif sumk(Qa,1) < b && b <= sumk(Qa,2) %BACKLOG WILL INCREASE 1 BECAUSE THERE IS COLLISION
                n=n+1; 
                packets_arriving(i)=aux+1;
                packets_leaving(i)=0;
            elseif sumk(Qa,2) < b && b <= sumk(Qa,3) %2 ARRIVALS, BACKLOG INCREASES TOO BECAUSE OF COLLISION
                n=n+2;
                packets_arriving(i)=aux+2;
                packets_leaving(i)=0;
            elseif sumk(Qa,3) < b && b <= sumk(Qa,4) %AND SO ON...
                n=n+3;
                packets_arriving(i)=aux+3;
                packets_leaving(i)=0;
            elseif sumk(Qa,4) < b && b <= sumk(Qa,5)
                n=n+4;
                packets_arriving(i)=aux+4;
                packets_leaving(i)=0;
            elseif sumk(Qa,5) < b && b <= sumk(Qa,6)
                n=n+5;
                packets_arriving(i)=aux+5;
                packets_leaving(i)=0;
            elseif sumk(Qa,6) < b && b <= sumk(Qa,7)
                n=n+6;
                packets_arriving(i)=aux+6;
                packets_leaving(i)=0;
            elseif sumk(Qa,7) < b && b <= sumk(Qa,8)
                n=n+7;
                packets_arriving(i)=aux+7;
                packets_leaving(i)=0;
            elseif sumk(Qa,8) < b && b <= sumk(Qa,9)
                n=n+8;
                packets_arriving(i)=aux+8;
                packets_leaving(i)=0;
            elseif sumk(Qa,9) < b && b <= sumk(Qa,10)
                n=n+9;
                packets_arriving(i)=aux+9;
                packets_leaving(i)=0;
            elseif sumk(Qa,10) < b && b <= 1
                n=n+10;
                packets_arriving(i)=aux+10;
                packets_leaving(i)=0;
            end
        end
    end
 backlog_array(i) = n; %Fill the array of the system's backlog 
end
 
figure(1) %Setting up the plotting environment for the backlog of the system
plot(slots,backlog_array)
xlabel('Slot number')
ylabel('Backlog of the system')
title('Backlog of the system vs. Slot number')
 
figure(2) %Setting up the plotting environment for the packets of the system
plot(slots,packets_arriving)
hold on
plot(slots,packets_leaving,'k')
grid on
xlabel('Slot number')
ylabel('Number of packets')
title('Number of packets entering (blue) or leaving (black) the system vs. Slot number')
 
figure(3) %Setting up the plotting environment for the histogram of the states
max = max(backlog_array);
hist(backlog_array,max) %array with the counts of the times in each state
nelements = hist(backlog_array,100); %In each bin, it counts how many times element i was seen ("numerical plot of the hist function")
for i = 1:m
    state_probs(i) = nelements(i)/1000;
end
xlabel('Number of state (n)')
ylabel('Frequency of each state')
title('Frequency of each state n')
 
%Calculation of the simulated probability of success:
ps_simulation = success/t;
%Calculation of the theoretical probability of success:
for i = 1:m
    aux = (binopdf(1,m-(i-1),qa)*binopdf(0,(i-1),qr))+(binopdf(0,m-(i-1),qa)*binopdf(1,i-1,qr));
    aux2 = aux * state_probs(i);
    ps_theory = ps_theory + aux2;
end
 
%Attempt rate plot (theoretical value):
for n = 1:m
    Att_Rate(n) = (m-(n-1))*qa + (n-1)*qr;
end
figure(4)
plot(states,Att_Rate,'r')
axis([0 100 0 1])
xlabel('State, n')
ylabel('Attemp rate, G(n)')
title('Attempt rate as a function of the current state n')
 
disp(['Simulated Probability of Success: ' num2str(ps_simulation) '.'])
disp(['Theoretical Probability of Success: ' num2str(ps_theory) '.'])