function [t,Ni,Ei] = chap_8_prob_14(tFinal,lmean,mmean)
% 
% This code simulates a M/M/1 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 given by an exponential distributed
% with rate mu
%
% 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
%
% 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  = 1/0.2;
end
lambda = 1/lmean;
if( nargin < 3 ) 
  % mu the rate of our exponential random variable: 
  mmean = 1/0.1;
end
mu = 1/mmean; 

% 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 )
  tia = (-1/lambda)*log( rand(1) );   % <- generate an interrival time  
  ti  = [ ti, tia ];                  % <- append it
  tsa = (-1/mu)*log( rand(1) );       % <- generate this customers service time
  ts  = [ ts, tsa ]; 
  tSum = tSum + tia;                  % <- compute the total time to have all these customers ARIVE
  nCusts = nCusts + 1 ;  
end
tcs = cumsum( 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 ];
    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 ];
      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 ); 

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

return;