function [AAvgMinutesInBank, BAvgMinutesInBank, AvgNumHoursServerGoesHome, RatioAToB] = BankQueue()
% Example: [AAvgMinutesInBank, BAvgMinutesInBank, AvgNumHoursServerGoesHome, RatioAToB] = BankQueue()

for runCount=1:100
  
    clear AWaitingTimes
    clear AInterArrivalTimes
    clear AServiceTimes
    clear ADepartures
    clear ATimeInBank
    clear AArrivalTimes
    clear ADepartureTimes
    
    clear BWaitingTimes
    clear BInterArrivalTimes
    clear BServiceTimes
    clear BDepartures
    clear BTimeInBank
    clear BArrivalTimes
    clear BDepartureTimes
    
    clear Queue
    
    % Customer Type A:
    i=1;
    AWaitingTimes(i)=0;  % Initial condition
    AInterArrivalTimes(i)=exprnd(6);
    AServiceTimes(i)=exprnd(3);
    ADepartures(i)=sum(AInterArrivalTimes(1:i))+AServiceTimes(i)+AWaitingTimes(i);
    ATimeInBank(i)=AServiceTimes(i);
    while sum(AInterArrivalTimes)<480  % 8 hours * 60 minutes/hour
        i=i+1;
        AInterArrivalTimes(i)=exprnd(6); % T(j)~EXP(1/lambda), A(j)=sum(T(1:j))
        AServiceTimes(i)=exprnd(3); % V(j)~EXP(lambda)
        AWaitingTimes(i)=max(0, AWaitingTimes(i-1)+AServiceTimes(i)-sum(AInterArrivalTimes(1:i)));
        ADepartures(i)=sum(AInterArrivalTimes(1:i))+AServiceTimes(i)+AWaitingTimes(i); % D(j)=sum(T(1:j))+V(j)+W(j)
        ATimeInBank(i)=AInterArrivalTimes(i)+AWaitingTimes(i)+AServiceTimes(i);
    end
    
    
    % Customer Type B:
    j=1;
    BWaitingTimes(j)=0;  % Initial condition
    BInterArrivalTimes(j)=exprnd(10);
    BServiceTimes(j)=exprnd(10);
    BDepartures(j)=sum(BInterArrivalTimes(1:j))+BServiceTimes(j)+BWaitingTimes(j);
    BTimeInBank(j)=BServiceTimes(j);
    while sum(BInterArrivalTimes)<480  % 8 hours * 60 minutes/hour
        j=j+1;
        BInterArrivalTimes(j)=exprnd(10); % t(j)~EXP(1/lambda), A(j)=sum(T(1:j))
        BServiceTimes(j)=exprnd(10); % V(j)~EXP(lambda)
        BWaitingTimes(j)=max(0, BWaitingTimes(j-1)+BServiceTimes(j)-sum(BInterArrivalTimes(1:j)));
        BDepartures(j)=sum(BInterArrivalTimes(1:j))+BServiceTimes(j)+BWaitingTimes(j); % D(j)=sum(T(1:j))+V(j)+W(j)
        BTimeInBank(j)=BInterArrivalTimes(j)+BWaitingTimes(j)+BServiceTimes(j);
    end
    
    
    % Concatenate customer A queue and customer B queue and create stairs-plot:
    AArrivalTimes=[1.*ones(1,(length(AInterArrivalTimes)-1)); cumsum(AInterArrivalTimes(1:i-1))];
    BArrivalTimes=[1.*ones(1,(length(BInterArrivalTimes)-1)); cumsum(BInterArrivalTimes(1:j-1))];
    ADepartureTimes=[-1.*ones(1,(length(ADepartures)-1)); ADepartures(1:i-1)];
    BDepartureTimes=[-1.*ones(1,(length(BDepartures)-1)); BDepartures(1:j-1)];
    Queue=sortrows([AArrivalTimes BArrivalTimes ADepartureTimes BDepartureTimes]',2);
    Queue(:,1)=cumsum(Queue(:,1));
    
    AMinutesInBank(runCount)=mean(ATimeInBank);
    BMinutesInBank(runCount)=mean(BTimeInBank);
    LastDepartureTime(runCount)=Queue(end,2)/60;
    TypeAToTypeB(runCount)=length(AArrivalTimes)/length(BArrivalTimes);
end

stairs(Queue(:,2),Queue(:,1))
axis([0 max(Queue(:,2))+1 -1 max(Queue(:,1))+1])
xlabel('Time (in minutes)')
ylabel('Number of Customers in the Queue')


AAvgMinutesInBank=mean(AMinutesInBank);
BAvgMinutesInBank=mean(BMinutesInBank);
AvgNumHoursServerGoesHome=mean(LastDepartureTime);
RatioAToB=mean(TypeAToTypeB);