% Parametry.
lambda = 4.5;    % intenzita příchodů (průměrný počet za hodinu)
mu = 5;          % intenzita obsluhy (průměrný počet za hodinu)
t_max = 500;     % délka simulace (hodin) - po vyzkoušení zvětšit

% Inicializace proměnných.
t = 0;           % aktuální simulační čas
fronta = 0;      % aktualní délka fronty
obsluha = 0;     % =1, pokus se právě někdo obsluhuje
t_obsluha = inf; % čas nejbližší obsluhy, inf = "infinity"
t_prichod = -1/lambda*log(rand(1)); % čas nejbližšího příchodu

% Inicializace proměnných, do kterých se zapisují simulační data.
prichody = NaN;  % časy odchodů, NaN = "not a number"
odchody = NaN;   % časy odchodů
casy = NaN;      % časy jednotlivých událostí
vsystemu = NaN;  % počty lidí v časech jednotlivých událostí

% Simulace.
while t < t_max
   % Tato cast slouzi jen pro kresleni grafu.
   casy(end+1) = t;
   vsystemu(end+1) = fronta + obsluha;
   
   % Tato cast obsahuje kod pro vlastni simulaci.
   if t_prichod < t_obsluha
      t = t_prichod;
      t_prichod = t - 1/lambda*log(rand(1));
      prichody(end+1) = t;
      if obsluha == 0
         obsluha = 1;
         t_obsluha = t - 1/mu*log(rand(1));
      else
         fronta = fronta + 1;
      end
   else
      t = t_obsluha;
      odchody(end+1) = t;
      if fronta == 0
         obsluha = 0;
         t_obsluha = inf;
      else
         fronta = fronta - 1;
         t_obsluha = t - 1/mu*log(rand(1));
      end
   end
end

% Textove vyhodnoceni vysledku.
obslouzeno = size(odchody,2) - 1;
clc
doby = odchody(2:obslouzeno+1) - prichody(2:obslouzeno+1);
disp(' ')
disp(['             Delka simulace: ' num2str(t_max)])
disp(['         Intenzita prichodu: ' num2str(lambda)])
disp(['          Intenzita odchodu: ' num2str(mu)])
disp(' ')
disp(['Celkem obslouzeno pozadavku: ' num2str(obslouzeno)])
disp(['V systemu na konci simulace: ' num2str(fronta+obsluha)])
disp(['    Prumerna doba v systemu: ' num2str(mean(doby))])

% Histogram rozdělení doby strávené v systému (podle požadavků). 
subplot 211
pocetint = 50; % počet intervalů, do kterých se seskupí lidé podle doby
x = max(doby)/pocetint:max(doby)/pocetint:max(doby); % meze intervalů
hist(doby,x) % vlastní graf
title('Rozdělení doby strávené v systému (podle požadavků)', ...
    'FontWeight','bold')
xlabel('doba v systému')
ylabel('počet lidí')

% Graf pro vývoj počtu lidí v systému.
subplot 212
stairs(casy(2:end), vsystemu(2:end))
title('Vývoj počtu lidí v systému','FontWeight','bold')   
xlabel('čas')
ylabel('aktuální počet lidí')