function[a,PZ,PQ,lambdaP,lambdaZ,kappa,ENS,ENQ,EN,EWS,EWQ,EW]=odmitani(lambda,mi,n,m);
% [a,PZ,PQ,lambdaP,lambdaZ,kappa,ENS,ENQ,EN,EWS,EWQ,EW]=odmitani(lambda,mi,n,m)
%
% Vypocita stacionarní rozlozeni, vyuziti a charakteristiky systemu 
% hromadne obsluhy s omezenou kapacitou M|M|n|m|FIFO s odmitanim.
%
% Vstupní parametry:
% lambda .... parametr vstupniho proudu
% mi ........ parametr obsluhy
% n ......... pocet linek obsluhy
% m ......... kapacita systému
%
% Vystupní parametry:
% a ......... stacionární rozlození
% PZ ........ pravdepodobnost, ze prichozi zakaznik bude odmitnut
% PQ ........ pravdepodobnost, ze prichozi zakaznik bude cekat ve fronte
% lambdaP ... stredni hodnota poctu prijatych zakazniku za jednotku casu
% lambdaZ ... stredni hodnota poctu odmitnutych zakazniku za jednotku casu
% ENS ....... stredni hodnota poctu obsluhovanych zakazniku
% ENQ ....... stredni hodnota poctu zakazniku ve fronte
% EN ........ stredni hodnota poctu zakazniku v systemu
% EWS ....... stredni hodnota doby, kterou zakaznik stravi obsluhou
% EWQ ....... stredni hodnota doby, kterou zakaznik stravi ve fronte
% EW ........ stredni hodnota doby, kterou zakaznik stravi v systemu

beta=lambda/mi;
ro=beta/n;
a(1)=inv(sum(beta.^(0:n-1)./factorial(0:n-1))+n^n/factorial(n)*sum(ro.^(n:m)));
a(2:n+1)=(beta.^(1:n)./factorial(1:n))*a(1);
a(n+2:m+1)=(n^n/factorial(n))*ro.^(n+1:m)*a(1);
PZ=a(m+1);
if ro==1
  PQ=a(n+1)*(m-n);
else
  PQ=a(n+1)*(1-ro^(m-n))/(1-ro);
end;
lambdaP=lambda*(1-PZ);
lambdaZ=lambda*PZ;
kappa=ro*(1-PZ);
ENS=beta*(1-PZ);
ENQ=sum(((n+1:m)-n).*a(n+2:m+1));
EN=ENQ+ENS;
EW=EN/lambda;
EWS=ENS/lambda;
EWQ=ENQ/lambda;