function [aveNwait, lastNwait, avwaitT] = system2(N,mc,ms,k,accu)
%SYSTEM2 simulates a checkout system with shortest queue allocation of customers.  
%Created for 1.010 Uncertainty in Engineering
%CC Andreas Langousis and Daniele Veneziano
%MIT Department of Civil and Environmental Engineering
%Additional code by Violeta Ivanova, violeta@mit.edu

%USAGE system2(N,mc,ms,k,accu)
%N:number of counters
%mc:mean customer arrival time 
%ms:mean service time at counters
%k:number of system states to be simulated
%accu:accuracy. The program stops when either k is reached or the 
%coefficient of variation is les than accu

%VARIABLES USED IN THE PROGRAM
%coall:line matrix of dimention N that includes the total number of
%customers that have been served and wait to be served at each counter
%cowait:line matrix of dimention N that includes the number of
%customrers that wait to be served at each counter
%cowaitave:line matrix of dimention N that includes the average number of
%customrers that wait to be served at each counter
%cowaitstdev:line matrix of dimention N that includes the stdev of the
%number of customrers that wait to be served at each counter
%waitT:matrix with N columns and variable lines that includes the waiting
%time of each customer at each counter
%n:number of counters that are working at the current state
%c:optimization parameter of the system (c>1 no infinite queues are formed)
%mt:mean time to system status change
%P:Probability of a new customer arrival at the current state
%t:random time to next system status change
%flag: 0 if a customer left, 1 if a new customer arrived
%g:counter that changes status
%avwaitT:matrix with 3 columns and variable # of lines. 1st column is the
%simulation step, 2nd column real time, 3rd column average waiting time.
%v:coefficient of vriation
%temp:temporary matrix needed for taking randomly the minimum between equal values.
%matlab routine for the min takes the first minimum on a line matrix.
%temp1:Temporary storage matrix

disp('System 2: Shortest Queue - upon arrival customers go to the shortest queue.'); 
disp('Simulations in progress ...');

%define minimum number of trials for the coefficient of variation check
if k<=1500;
    k=1500;
end

coall=zeros(1,N);
cowait=zeros(1,N);
cowaitave=zeros(1,N);
cowaitstdev=zeros(1,N);
temp=zeros(1,N);
waitT=zeros(k,N);
avwaitT=zeros(k+1,3);
check=zeros(1,1);
temp1=zeros(1,1);
n=0;
i=0;

%estimate optimization parameter c
c=N*mc/ms;
v=1000;

%starting simulation of the system
while(and((i<k),(v>=accu)));
    i=i+1;
    %estimate mean time to system status change
        mt=1/(1/mc+n/ms);
        P=1/(1+n*mc/ms);
        %generate random numbers and define the state of the system
        t=exprnd(mt);
        if rand(1)<= P;
            flag=1;
        else
            flag=0;
        end
            
        if flag==1;
            %find counter that gained a customer
            %random choice of the minimum
            for l=1:N
                temp(l)=cowait(l)+rand(1)*0.001;
            end
            [e,g]=min(temp);
            if cowait(g)==0;
                n=n+1;
            end
            cowait(g)=cowait(g)+1;
            coall(g)=coall(g)+1;
            
            %increase waiting times
            for l=1:N;
                if l==g 
                    for h=1:cowait(l)-1;
                        waitT(coall(l)-h,l)=waitT(coall(l)-h,l)+t;
                    end
                else
                    for h=0:cowait(l)-1;
                        waitT(coall(l)-h,l)=waitT(coall(l)-h,l)+t;
                    end
                end   
            end
        else
            %find counter that lost a customer
            g=0;
            while(g==0);
                g=round(rand(1)*N);
                if g==0;
                  g=N;
                end
                if cowait(g)==0;
                    g=0;
                end
            end
            
            %increase waiting times
            for l=1:N;
                for h=0:cowait(l)-1;
                    waitT(coall(l)-h,l)=waitT(coall(l)-h,l)+t;
                end
            end
            
            cowait(g)=cowait(g)-1;
            if cowait(g)==0;
                n=n-1;
            end
        end
        
        %estimating average waiting time
        avwaitT(i+1,1)=i;
        avwaitT(i+1,2)=avwaitT(i,2)+t;
        
        s=0;
        for l=1:N;
            s=s+coall(l)-cowait(l);
        end
        for l=1:N;
            for h=1:coall(l)-cowait(l);
               avwaitT(i+1,3)=avwaitT(i+1,3)+waitT(h,l)/s;
            end
        end
        
        %variation coefficient
        check(i+1)=avwaitT(i+1,3);
        if i>1000
            v=((var(check))^0.5)/mean(check);
        end
        for j=1:N 
           temp1(i,j)=cowait(j);
        end
           
end
for j=1:N
    cowaitave(j)=mean(temp1(:,j));
    cowaitstdev(j)=var(temp1(:,j))^0.5;
end

%column vector of average number of customers in front of counters
aveNwait = cowaitave';
disp('Average number of customers at the counters: ');
disp(cowaitave);

%column vector of the stdev of number of customers at each counter 
lastNwait = cowaitstdev';

disp('Standard deviation of the number of customers at the counters:');
disp(cowaitstdev);

system2=avwaitT(i+1,3);

return