%--------------------------------------------------------------------------
% Simulation of an M/M/1 queue
% Packets' inter-arrival time distribution is poisson with parameter lambda
% Packets' service time distribution is exponential with mu=1
% Output: average delay and queue length
%-------------------------------------------------------------------------
clc;
clear all;
close all;

npkt= 20000;    % number of pkts which are served
mu= 1.0;        % Service rate


lambda= 0.01:0.02:0.96; % rho=0:0.05:0.96

% Track the system state and pkt delay for each rho
N= zeros(length(lambda),1);
T= zeros(length(lambda),1);

for k=1:length(lambda)
    
    now= 0;                          % current time
    systemState= zeros(npkt,1);      % number of pkts in the system at each arrival
    departureTime= zeros(npkt,1);    % departure time of each pkt
    npktServed= 0;
    
    % Generate pkts inter-arrival and service distribution time
    rng('shuffle');
    interarrivalTime= exprnd(1/lambda(k),[npkt 1]);
    rng('shuffle');
    serviceTime= exprnd(1/mu,[npkt 1]);
    
    % the first pkt arrival and initial next departure
    nextdeparture= inf;
    nextarrival= interarrivalTime(1);
    npktArrived= 1;
    currentLength= 0;
    
    while npktServed < npkt
        if nextarrival < nextdeparture
            now= nextarrival;
            systemState(npktArrived)= currentLength;
            currentLength= currentLength + 1; % increment number of pkts in the system
            if currentLength == 1,
                nextdeparture= now + serviceTime(npktServed+1);
            end
            if npktArrived < npkt
                npktArrived= npktArrived+1;
                nextarrival= now + interarrivalTime(npktArrived); % update next arrival
            else
                nextarrival = inf;
            end
        else
            now= nextdeparture;
            currentLength= currentLength - 1;
            npktServed = npktServed+1;
            departureTime(npktServed)= now;
            if currentLength > 0,
                nextdeparture= now + serviceTime(npktServed+1);
            else
                nextdeparture= inf;
            end
            
        end
    end
    
    N(k)= mean(systemState);
    T(k)= mean(departureTime - cumsum(interarrivalTime));
    
end

% plot analytical and simulation mean delay
h1=figure(1);
set(h1,'color',[ 1 1 1]);
hold on;
Tanalysis = 1./(mu-lambda);
plot(lambda,Tanalysis,'r--','LineWidth',2);
plot(lambda,T,'k','LineWidth',2);
xlabel('\rho');
ylabel('Delay');
legend('Analysis','Simulation');

% plot analytical and simulation number of pkts in the system
h2=figure(2);
set(h2,'color',[ 1 1 1]);
hold on;
Nanalysis = lambda./(mu-lambda);
plot(lambda,Nanalysis,'r--','LineWidth',2);
plot(lambda,N,'k','LineWidth',2);
xlabel('\rho');
ylabel('Nubmer of pkts in system');
legend('Analysis','Simulation');

% output process of the queue
h3=figure(3);
set(h3,'color',[ 1 1 1]);
cdfplot(diff(departureTime));
hold on;
cdfplot(interarrivalTime);