function [d,u,nq] = queue2(n,k,lambda,mu)
% Simulate delays in queue for queuing model (M/M/k)
% ------------------------------------ 
% n the number of customers
% k the number of servers
% lambda the arrival rate of customers
% mu the service rate
% ------------------------------------
% d the delay times for each customer
% u the utilization of servers
% nq the average queue length
% ------------------------------------

t = 0;                % Simulation clock
ta = exprnd(lambda);  % The time of next customer arrival
td = inf;             % Departure time of customers (set to a dummy value)

served = [];          % Customers being served
queue = [];           % Queued customers

n_in = 0;             % Number of customer arrivals
n_out = 0;            % Number of customer departures

d = [];               % Service delays
u = 0;                % Utilization of server
                      % (time interval * fraction of servers having
                      % customer in that time interval)
nq = 0;               % Number in queue statistic
                      % (time interval * length of queue)

% The main simulation loop
while n_out < n
    
    % Next event is arrival
    if ta < td                    
        
        % Update statistics
        
        % Number of customers entered the system
        n_in = n_in + 1;    
        % Utilization statistic
        u = u + (ta-t) * length(served) / k;
        % Number of customers in queue statistic
        nq = nq + (ta-t) * length(queue);

        % Add either to server or queue
        if length(served) < k
            % Generate departure time, add to served vector and sort
            % so that earliest departure is the first element
            served = sort([served ta+exprnd(mu)]);
            % Customer went straight to server => no delay
            d(end+1)= 0;
        else 
            % Add the arrival time to the queue vector 
            queue = [queue ta];
        end
        
        % Update simulation clock 
        t = ta;
        % Generate next arrival time
        ta = ta + exprnd(lambda);
        % Set next departure time as the earliest departure time, i.e.,
        % the first element of the served vector
        td = served(1);
        
        
    else
        % Next event is departure
        
        % Update statistics
        n_out = n_out + 1;
        u = u + (td-t) * length(served) / k;
        nq = nq + (td-t) * length(queue);
        
        % Departure, remove first element of served vector
        served = served(2:end);

        % Take next customer from queue, if queue is not empty
        if ~isempty(queue)
            % Add new customer to served list and generate departure time
            served = sort([served td+exprnd(mu)]);
            % Record delay: Now - time the customer entered the queue
            d(end+1)= td - queue(1);
            % Remove the customer, i.e., the first element from the queue
            queue = queue(2:end);
        end
        
        % Update simulation clock
        t = td;
        % If servers are empty
        if isempty(served)
            % Set next departure time high so that the event cannot happen
            td = inf;
        else
            % If there is customers in the servers, set the next departure
            % time as the earliest departure time in served
            td = served(1);
        end
        
    end
end

% Calculate utilization: In u, we have the added time intervals
% between the events times the fraction of servers having customer in the
% intervals. To get utilization we simply divide u with the simulation time 
% in the end of the simulation
u = u / t;
% Calculate average queue length: Similar to utilization
nq = nq / t;