% Simple queuing theory simulation, M/M/1 queue
% Single server, single queue

a = 1;  % average number of arrivals per minute
b = 1.5;  % average number of people served per minute

ncust = 1000;
% Notation:
%  at = arrival time of a person joining the queue
%  st = service time once they reach the front
%  ft = finish time after waiting and being served.
%

% initialize arrays:
at = zeros(ncust,1);  
ft = zeros(ncust,1);

% Generate random arrival times assuming Poisson process:
r = rand(ncust,1);
iat = -1/a * log(r);   % Generate inter-arrival times, exponential distribution
at(1) = iat(1);        % Arrival time of first customer

for i=2:ncust
   at(i) = at(i-1) + iat(i);   % arrival times of other customers
   end

% Generate random service times for each customer:
r = rand(ncust,1);
st = -1/b * log(r);   % service times for each


% Compute time each customer finishes:
ft(1) = at(1)+st(1);  % finish time for first customer
for i=2:ncust
   % compute finish time for each customer as the larger of
   %   arrival time plus service time (if no wait)
   %   finish time of previous customer plus service time (if wait)
   ft(i) = max(at(i)+st(i),  ft(i-1)+st(i));
   end

total_time = ft - at;           % total time spent by each customer
queue_time  = total_time - st;   % time spent waiting before being served

ave_service_time = sum(st)/ncust
ave_queue_time = sum(queue_time)/ncust
ave_total_time = sum(total_time)/ncust

% compute length of queue as a function of time.
% merge together the arrival and finish times and sort them
% into one list, keeping track of whether queue length changes by +1
% or -1 at each time:

event_times = [at; ft];

queue_change = [ones(ncust,1); -ones(ncust,1)];
[sorted_times, i_times] = sort(event_times);
sorted_change = queue_change(i_times);
queue_length = zeros(2*ncust,1);
for i=2:(2*ncust)
   queue_length(i) = queue_length(i-1) + sorted_change(i-1);
   end

figure(1)
plot(sorted_times,queue_length)
h1 = title('queue length vs. time');
set(h1,'fontsize',15)


% Plot histogram of waiting times:

figure(2)
hist(total_time,0:.5:20)
h1 = title('histogram of waiting times');
set(h1,'fontsize',15)
axis([-.5 20.5 0 400])
h1=text(6,380,['a =',num2str(a)]);
set(h1,'fontsize',15)
h1=text(6,350,['b =',num2str(b)]);
set(h1,'fontsize',15)
h1=text(6,320,['average total time =',num2str(ave_total_time)]);
set(h1,'fontsize',15)
h1=text(6,290,['average service time =',num2str(ave_service_time)]);
set(h1,'fontsize',15)