%
%  macm 202 -- 11 feb 03 -- djm
%
%  w6queue.m:  simulate queueing by exponential transition process
%
%     - use direct exponential distribution
%     - multiple servers
%     - single-queue or multi-queues!
clear

study = 1

%  simulation parameters & pseudorandom seed
Nevents = 1000;  Nservers = 3;  queue_max = 30;
muQ = 2.5;  muS = 1/1;  big = 10^6;

rand('state',7654321);

%  initialize data structures
event_data    =     zeros(Nevents,7);

queue_size    =     zeros(Nservers,1);
queue_array   = -big*ones(Nservers,queue_max);
server_array  =     zeros(Nservers,1);
server_status =      ones(Nservers,1);

clock         = 0;
next_arrival  = -log(rand)/muQ;

%  simulation loop
for k=1:Nevents;
    [interval eventNo] = ...
        min([server_array+big*server_status ; next_arrival]);
    
    %  update all times
    clock             = clock + interval;
    event_data(k,1:2) = [eventNo clock];

    queue_array       = queue_array + interval;

    server_array      = server_array - interval;
    next_arrival      = -log(rand)/muQ;

    %  process event
    switch eventNo
        %  if queue arrival
        case{Nservers+1}
            %  check if sleeping server
            [work,serverNo] = max(server_status);
            if(work==1)
                %  go straight to the wicket
                server_array(serverNo)  = -log(rand)/muS;
                server_status(serverNo) = 0;
                event_data(k,6)         = server_array(serverNo);
            else
                %  gotta wait (one line, or shortest line!)
                switch study
                    case{1}
                        shortNo = 1;
                    case{2}
                        [junk,shortNo] = min(queue_size);
                end
                queue_size(shortNo) = queue_size(shortNo) + 1;
                queue_array(shortNo,queue_size(shortNo)) = 0;
                event_data(k,3:4)   = [shortNo queue_size(shortNo)];
            end

        %  if service completion    
        otherwise
            %  check for anyone waiting
            if (sum(queue_size)>0)
                %  if your line empty, most impatient from other queues!
                if (queue_size(eventNo)==0)
                    [junk,nextNo] = max(queue_array(:,1));
                else
                    nextNo = eventNo;
                end
                %  reset server task & shift queue
                server_array(eventNo) = -log(rand)/muS;           
                event_data(k,6)       = server_array(eventNo);

                queue_size(nextNo)    = queue_size(nextNo) - 1;
                event_data(k,5)       = queue_array(nextNo,1);

                queue_array(nextNo,1:queue_size(nextNo)) = ...
                    queue_array(nextNo,2:queue_size(nextNo)+1);
                queue_array(nextNo,queue_size(nextNo)+1) = -big;

                event_data(k,3:4)     = [nextNo queue_size(nextNo)];
            else
                %  nobody in queue
                server_array(eventNo)  = 0;
                server_status(eventNo) = 1;
            end
    end
    event_data(k,7) = sum(queue_size);
end    

%  plot queue_size vs time
figure(study);  clf
time  = [[0 ; event_data(1:end-1,2)] event_data(1:end,2)]';  
time  = time(:);
queue = [event_data(1:end,7) event_data(1:end,7)]';  
queue = queue(:);

plot(time,queue)
title('\bf queue size vs time')
xlabel('\bf clock time')
ylabel('\bf queue size')