function computeQueueCapacity(phase, Kcheck)
% finds queue capacity K so that the blocking rate is smaller than the 
% given value Bmax
Bmax = 1e-3;
T = 1e5;               % simulation time for root finding
Tmax = 1e7;            % simulation time for checking
K0 = 15;               % too small capacity
K1 = 35;               % too large capacity

if phase == 1
  % find reasonable values
  % auxiliary function to solve for B = Bmax
  funk = @(K) blockingRate(round(K), T) - Bmax;
  fzero(funk, [K0, K1])   
  % stop manually when result fluctuates -> K = 27/28
elseif phase == 2
  % check with large T
  K = Kcheck;
  blockingRate(K, Tmax);
  blockingRate(K+1, Tmax);
  blockingRate(K+2, Tmax);
end

function B = blockingRate(K, T)
% gets blocking rate B for given queue capacity K by running the model

% list all dialog parameters with
%    get_param('queuecapacity1C/Queue','DialogParameters')
set_param('queuecapacity1C/Queue','Capacity', num2str(K))
simOut = sim('queuecapacity1C', 'StopTime', num2str(T));
B = simOut.nBlockings(end)/T;
fprintf('run: K = %3d,  B = %8.6f\n', K, B)
