%% Simulace telefonní linky.
lambda = 1;     % Intenzita příchodů hovorů.
mu = 3;         % Intenzita odbavení hovorů.

hovory = 1000;  % Počet hovorů (volíme).
volani = 0;     % Počet volání (zjistíme, zatím 0).
t_prichod = 0;  % Čas nejbližšího příchodu.
t = 0;          % Aktuální simulační čas.

for i = 1:hovory
    t = t - log(rand)/mu;    % Čas ukončení telefonátu.
    while t_prichod <= t     % Je-li linka obsazena, 
        volani = volani + 1; % zavěsím.
        t_prichod = t_prichod - log(rand)/lambda;  
    end
    t = t_prichod;  % Čas začátku nového telefonátu.
end

volani
podil_zavesenych = (volani-hovory)/volani

%% Grafy pro telefonní linku (z přednášky).
Q = [-lambda lambda
          mu    -mu];
      
% Časová škála.
tcount = 100;
step = 0.02;
tmax = step*tcount;

% Výpočet P(t).
P = zeros(2,2,tcount+1);
P(:,:,1) = eye(2);
for t = 1:tcount
    P(:,:,t+1) = expm(t*step*Q);
end

% Uloží p_ij(t) do vektoru.
p00 = squeeze(P(1,1,:));
p01 = squeeze(P(1,2,:));
p10 = squeeze(P(2,1,:));
p11 = squeeze(P(2,2,:));
t = (0:step:tmax)';

% Grafy
plot(t,p00,'b',t,p01,'r',t,p10,'b:',t,p11,'r:')
ylabel('{\itp_{ij}}({\itt})','FontSize',12,'Rotation',0)
xlabel('\itt','FontSize',12)
legend('{\itp}_{00}({\itt})','{\itp}_{01}({\itt})','{\itp}_{10}({\itt})','{\itp}_{11}({\itt})')
title('\lambda = 1, \mu = 3')