function [aveNwait, lastNwait, avwaitT] = system3(N,mc,ms,k,accu)
%SYSTEM1 simulates a checkout system with a single 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 system3(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:total number of customers that have been served and wait to be served (refers to all
%counters)
%cowait: number of customers that wait to be served (refers to all
%counters)
%cowaitave:this is the average number of customers in the line
%cowaitstdev:this is the standard deviation of customers in the line
%waitT:single column matrix with variable lines that includes the waiting
%time of each customer
%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.
%temp:temporary storage matrix

disp('System 3: Single Queue - customers go to the first available counter.'); 
disp('Simulations in progress ...');

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

coall=0;
cowait=0;
cowaitave=0;
cowaitstdev=0;
waitT=zeros(k,1);
avwaitT=zeros(k,3);
check=zeros(1,1);
temp=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;
            %customer arrival
            if cowait<N;
                n=n+1;
            end
            cowait=cowait+1;
            coall=coall+1;
            
             %increase waiting times 
             for h=1:cowait-1;
                 waitT(coall-h)=waitT(coall-h)+t;
             end       
        else
            %customer departure
            %increase waiting times
            for h=0:cowait-1;
                 waitT(coall-h)=waitT(coall-h)+t;
            end
            cowait=cowait-1;
            if cowait<N;
                n=n-1;
            end
        end
       %estimating average waiting time
        avwaitT(i+1,1)=i;
        avwaitT(i+1,2)=avwaitT(i,2)+t;
        for h=1:coall-cowait;
           avwaitT(i+1,3)=avwaitT(i+1,3)+waitT(h)/(coall-cowait);
        end  
        %variation coefficient
        check(i+1)=avwaitT(i+1,3);
        if i>1000
            v=((var(check))^0.5)/mean(check);
        end
        temp(i)=cowait;
end
%print the average number of customers in the line
cowaitave=mean(temp);
cowaitstdev=var(temp)^0.5;

%column vector of average number of customers in the line
aveNwait = zeros(N, 1);
aveNwait(1, 1) = cowaitave;
disp('Average number of waiting customers : ');
disp(cowaitave);

%column vector of the stdev of number of customers in the line
lastNwait = zeros(N, 1);
lastNwait(1, 1) = cowaitstdev;
disp('Standard deviation of the number of customers waiting: ');
disp(cowaitstdev);

system3=avwaitT(i+1,3);

return