%
%  macm 202 -- 02 mar 03 -- djm
%
%  w9q.m:  single server queue
%
%     - use direct exponential distribution
%     - for single-queue waiting time
clear

study = 1;

%  simulation parameters & pseudorandom seed

Nevents = 60;  Nservers = 1;  queue_max = 30;  rand('state',1234567);
%Nevents = 10000;  Nservers = 1;  queue_max = 50;  rand('state',1234567);
muQ = 2.5;  muS = 3/1;  big = 10^6;

%  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      = next_arrival - interval;

    %  process event
    switch eventNo
        %  if queue arrival
        case{Nservers+1}
            next_arrival      = -log(rand)/muQ;
            %  check if sleeping server
            [work,serverNo] = max(server_status);
            if(work==1)
                %  go straight to the wicket
                server_array(serverNo) = -log(rand)/muS;
                event_data(k,6) = server_array(serverNo);
                server_status(serverNo) = 0;
            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(1);  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')

%  extract waiting times (only a wait if a service is recorded in data)

waits   = event_data(:,5);
service = event_data(:,6);
[service Iservice] = sort(service);
waits = waits(Iservice);
[junk,begin] = max(service>0);

wsort = sort(waits(begin:end));
[junk,begin] = max(wsort>0);