function [L, W]=ggm_queue(a1, a2, b1, b2, s, T)
%
% Monte Carlo simulation of G/G/m queue
% Interarrival times ~ U(a1, a2)
% Service time ~ U(b1, b2)
% s=number of servers
% T=length of simulation run
% The output consists of the following data (for a single run)
% L=Number of customers in system at time T
% W=Average sojourn time of customers who have arrived and departed at time T
%
% Variables used 
% t=clock 
% ta=time of next arrival 
% td(k) = time of next departure from server k, k=1..s
% n=number of customers in the system 
% custserved = vector(s) , shows which customer is currently served in server k, k=1..s
% first_in_queue = which customer is first in queue (if queue is empty then first_in_queue = 0)
% Customers are ordered in order of arrival (FCFS discipline)
% NA=number of arrivals so far
% ND=number of departures so far
% custdata = matrix, row i contains the arrival and departure times of customer i
%
             	 td=zeros(s,1);
             	 custserved=zeros(s,1);
             	 t=0; n=0; first_in_queue=0;
          	 ta=a1+rand(1)*(a2-a1);
         	 for k=1:s td(k)=1e12; end;
         	 NA=0; ND=0;
                    	 custdata=zeros(1,2); 
	while t <=T
    		[mintd, kmin]=min(td);
    		[maxtd, kmax]=max(td);
    		if ta<=mintd  	% next event is arrival
        			t=ta;
        			NA=NA+1;
        			custdata(NA,1)=t; 
        			n=n+1;
        			ta=t + a1+rand(1)*(a2-a1); %generate next arrival time
%
% If there are free servers (n<s) then maxtd = 1e12 because at least one server does not have a valid departure time. In this case kmax is the first server with td(k) =1e12 and the customer is moved to that server
        			if n<=s 
        				custserved(kmax)=NA;
        				td(kmax)=t+b1+rand(1)*(b2-b1); 
        			elseif n==s+1 
        				first_in_queue=NA; %if all servers are busy and queue is empty coustomer NA takes the first place in queue
        			end;
		else		% next event is departure from server kmin
        			t=mintd;
        			ND=ND+1;
%custserved(kmin) shows which customer is the one who is departing. The departure time of this customer is updated in table custdata.
                			custdata(custserved(kmin),2)=t;      		
             			n=n-1;
% If first_in_queue >0 then the first in queue goes into the server that just got free. If first_in_queue =0 then the server becomes idle (custserved=0 and td=1e12).
        			if first_in_queue==0
        				td(kmin)=1e12; 
        				custserved(kmin)=0;
        			else 
        				custserved(kmin)=first_in_queue;
        			     	td(kmin)=t+b1+rand(1)*(b2-b1); 
% If queue becomes empty after the first customer moves to the free server then first_in_queue=0 otherwise first_in_queue is increased by one (because the customer waiting after this customer is necessarily the one who came next).
				if n==s 
					first_in_queue=0;
				else
					first_in_queue=first_in_queue+1;
				end;
        			end;
        		end;
        		% Uncomment the next line to track the evolution of the system at each event for verification purposes
        		%t, n, NA, ND, ta, td, kmin, kmax, first_in_queue, custserved,custdata
	end;
	L=n;
	SW=0;
	for i=1:NA
		if custdata(i,2)>0 SW=SW+custdata(i,2)-custdata(i,1); end;
	end;
	W=SW/ND;
endfunction;