function [d,busy_time,total_time,ninq,idle_time,nidle] = computercenter(n,warmup,k,lambda,mu)
% Simulate a computer center with k computers, i.e., a (M/M/k)-queue

%Input
%   n          the number of jobs that are simulated
%   warmup     the lenght of warmup period (#jobs)
%   lambda     the mean inter-arrival time of jobs
%   mu         the mean completion time of jobs

%Output
%   d          the queuing delays of n jobs
%   busy_time  the share of the total time all computers are busy
%   total_time the average processing time of a job
%   ninq       the average number of jobs in queue
%   idle_time  the share of the total time all computers are idle
%   nidle      the average number of idle computers


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

warmup_t=0;         % Length of the warmup in simulated time

served=[];          % Jobs being completed (holds departure times)
queue=[];           % Queued jobs (holds arrival times)

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

d=[];               % Queueing delays of jobs in the system
busy_time=0;        % a) The fraction of time all computers are busy
total_time=0;       % b) Mean total time of jobs in the system
ninq=0;             % c) Average number of jobs in queue
idle_time=0;        % d) The fraction of time all computers are idle
nidle=0;            % e) Average number of idle computers


% The main simulation loop
while n_out < n
    
    % If next arrival time is less than next departure, 
    % -> next event is arrival
    if ta < td                  
        
        % Update the number of jobs that have entered the center
        n_in = n_in + 1;    

        % Record performance if warmup has ended
        if n_in > warmup
            
            % Job is not executed immediately on submission
            if length(served) == k
                busy_time = busy_time + (ta-t);
            end

            % Total time the jobs spend in queue or server
            total_time = total_time + (ta-t) * (length(queue) + length(served));  
            
            % Number of jobs awaiting execution.
            ninq = ninq + (ta-t) * length(queue);
            
            % Computer center is idle.
            if isempty(served)
                idle_time = idle_time + (ta-t);          
            end

            % Number of idle computers.
            nidle = nidle + (ta-t) * (k-length(served));   
            
        
        end
            
        % Add either to server or queue
        if length(served) < k
            % Calculate departure time and add to served list
            served = sort([served ta + exprnd(mu)]);
            % Queuing delay is zero
            d(end+1)= 0;
        else
            % Add job arrival time to queue list
            queue = [queue ta];
        end
        
        % Update simulation clock 
        t = ta;
        % Next job arrival time
        ta = ta + exprnd(lambda);
        % Next job departure time
        td = served(1);
        
        
    else 
        % Next event is departure
        
        % Update statistics (A job leaves the center)
        n_out = n_out + 1;
        
        % Record if warmup has ended
        if n_in > warmup
            
            % a)
            if length(served) == k
                busy_time = busy_time + (td-t);          
            end

            % b)
            total_time = total_time + (td-t) * (length(queue)+length(served));    

            % c)
            ninq = ninq + (td-t) * length(queue);           
            
            % d)
            if isempty(served)
                idle_time = idle_time+(td-t);           
            end

            % e)
            nidle = nidle + ( td-t) * (k-length(served));    
        
        end
        
        % Departure (Earliest departure time is removed from the list)
        served = served(2:end);

        % Take next customer from queue, if queue is not empty
        if ~isempty(queue)
            % Determine departure time for the new job and add to served
            served = sort([served td + exprnd(mu)]);
            % Delay is time spent in the queue
            d(end+1)= td - queue(1);
            % Remove from the queue
            queue = queue(2:end);
        end
        
        % Update simulation clock
        t = td;
        if isempty(served)
            % If center is empty, reset next departure time
            td = inf;
        else
            % if center is not empty, set the next departure time
            td = served(1);
        end
        
    end

    % Record length of warmup in time
    if n_in > warmup && warmup_t == 0
        warmup_t = t;
    end
    
end

% Return exactly n delays
d = d(1:n);

% Calculate the performance metrics and take into consideration the warmup 
% time
busy_time = busy_time/(t-warmup_t);
total_time = total_time/(n_in-warmup);
ninq = ninq/(t-warmup_t);
idle_time = idle_time/(t-warmup_t);
nidle = nidle/(t-warmup_t);