function [t,Ni,Ei,Ai,Di] = chap_8_prob_13(tFinal,lmean,ul,ur)
% 
% This code simulates a single-server queueing system up to a time t.  
% That is it generates a single sample path of the number of customers 
% in the system as a function of time up to time "t".  
%
% The cusomers arrive according to a Poisson process with rate lambda.
%
% The service time for each customer is uniformly distributed between a range
% [ul,ur]
%
% Outputs:
%   t:  the times when the number of people in the system changes
%   Ni: the number of people in the system at each of the times in the array t
%   Ei: an indicator +1 => customer arrival; -1 => customer departure
%   Ai: the TOTAL number of arrivals  by time t 
%   Di: the TOTAL number of departure by time t
%
% Written by:
% -- 
% John L. Weatherwax                2006-09-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the% above email
% address.
% 
%-----

% set the random seed: 
%
%randn('seed',0); rand('seed',0); 

% use default values: 
% 
if( nargin < 2 )
  % lambda = rate (1/mean) of our Poisson process: 
  lmean  = .3;
end
lambda = 1/lmean;
if( nargin < 3 ) 
  % the limits of the uniform distribution of service times: 
  ul=1; ur=3; 
end

% generate the customer INTERarrival times (enough to get to time "t"): 
% -- 
tSum = 0.0; 
nCusts = 0; 
ti = []; ts = []; % storage for interarrival/service times
while( tSum <= tFinal )
  % generate an interrival time: 
  %
  if( length(ti)<4 )
    % special case code to use the GIVEN uniform random variables:
    book_interrival_unif = [0.843, 0.157, 0.9, 0.912]; 
    t_indx               = length(ti)+1;
    tia                  = (-1/lambda)*log( book_interrival_unif(t_indx) );
  else
    % the normal default case: 
    tia = (-1/lambda)*log( rand(1) );   
    tia = 100;                        % <- if needed specify a time so great in the future that we have no arrivals ... 
  end
  ti  = [ ti, tia ];                  % <- append it
  
  % generate a service time: 
  % 
  if( length(ts)<4 )
    % special case code to use the GIVEN uniform random variables:
    book_service_unif = [0.774, 0.244, 0.817, 0.799];
    t_indx            = length(ts)+1;
    tsa               = ul + (ur-ul)*book_service_unif( t_indx ); 
  else
    % the normal default case:
    tsa = ul + (ur-ul)*rand(1);
  end

  ts  = [ ts, tsa ]; 
  tSum = tSum + tia;                   % <- compute the total time to have all these customers ARIVE
  nCusts = nCusts + 1 ;  
end
tcs = cumsum( ti ); 


ts, ti

t  = [ 0 ]; Ei = [ 0 ]; 
ca = 2; cs = 2; 
nxtArrTime   = ti(1);
nxtSerTime   = ti(1)+ts(1); 
processedAll = 0;
while( ~processedAll )
  if( nxtArrTime <= nxtSerTime )
    t  = [ t, nxtArrTime ];
    Ei = [ Ei, +1 ];                      % <- a ARRIVAL is indicated by +1 
    if( ca<=nCusts )
      nxtArrTime = tcs(ca); ca=ca+1;
    else
      nxtArrTime = +Inf;
    end
  else
    % there must be a waiting customer to be serviced ... if not revert to the old nxtSerTime
    if( sum(Ei)<0 )
      nxtSerTime = tcs(cs-1)+ts(cs-1);
    else
      t  = [ t, nxtSerTime ];
      Ei = [ Ei, -1 ];                    % <- a DEPARTURE is indicated by -1 
      if( cs<=nCusts )
        nxtSerTime = nxtSerTime + ts(cs); cs=cs+1;
      else
        nxtSerTime = +Inf;
      end
    end
  end
  processedAll = ((ca>nCusts) && (cs>nCusts));
end
Ni = cumsum( Ei ); 
Ai = Ei; Ai( find(Ai<0) )=0; Ai = cumsum(Ai);
Di = Ei; Di( find(Di>0) )=0; Di = abs(cumsum(Di)); 

% trim our arrays so that that end at time "tFinal": 
if( 1 ) 
  inds = find( t >= tFinal );
  t(inds) = []; Ni(inds) = []; Ei(inds) = []; Ai(inds) = []; Di(inds) = []; 
  t(end+1) = tFinal; Ni(end+1) = Ni(end); Ai(end+1) = Ai(end); Di(end+1) = Di(end); 
end

return;