AAaaaaAAaaaaaAAa

function SP_LR_4_KHMELEV
  % Variant 16

  lambda = 8.208;
  mu = 3.8;

  t = []; % moment vremeni
  t_min(1) = t_obsl(1) = exprnd(1/mu); % vremya obslujivania zayavki
  t_ojz(1) = exprnd(1/lambda); % vremya ojidania sleduyuschey zayavki
  k = 0; % chislo zanyatyh priborov ravno kolichestvu zayavok, nahodyaschihsya v sisteme
  Type(1) = 1; % tip sobytia: '1' - poyavlenie novoy zayavki (C(L) += 1, apps(free) = 1); '2' - zavershenie obslujivania zayavki (C(L) -= 1, apps(busy) = 0)
  C(1) = 1; % sostoyanie sistemy
  free = 0;
  busy = 1;
  j(1) = J = 1;
  freed = J;
  doobsl = processed = [];

  % task 4:
  po = lambda/mu; % teoreticheskoe sredneye chislo zanyatyh priborov
  r_and_k0 = stationDist(po);
    r = r_and_k0(1:end-2) % da, zdes' minus 2
    k0 = r_and_k0(end) - 2; % i tut toje, dolgo ob"yasnyat'
    r(k0+1); % vsego u nas mojet byt' (k0+1) zanyatyh priborov, tak kak chitaem ne s 0, a s 1
  %print_it(r);
  apps = zeros(1, k0+1); % v nachale u nas net ni odnogo zanyatogo pribora



  %task 5:
  DELTAS = [];
  t(1, 1) =  t_z(1) = t_sob(1) = exprnd(1/lambda);

  t_kob(1) = t_z(1) + t_obsl(1);
  apps(1) = 1; % zanyali pervy pribor
  L = 1; % nomer sobytia
  k(1) = 1;

  c = C(L); m = 0;
  delta = max( t(:, 1)' ./ t_sob(1) .- r );
  t_sob(L) = t_ojz(1);
  while delta >= 0.00001 && c < k0 && m < 100
    m += 1;
    delta_or_C = modeling;
    delta = delta_or_C(1);
    c = delta_or_C(2);

    DELTAS(end+1) = delta;
  endwhile
  K = length(DELTAS);
  j



  %______________________________________used functions:
  function delta_or_C = modeling
  L += 1; % nomer sobytia L

  Type(L) = type(t_min(L-1), t_ojz(L-1)); % tip sobytia L
  printf('\n---\nt_ojz = %i\n', t_ojz(L-1));
  printf('t_min = %i\n', t_min(L-1));
  printf('tip = %i\n\n', Type(L));
  switch Type(L)
    case 1 % tip 1:
      t_sob(L) = t_sob(L-1) + t_ojz(L-1); % moment nastuplenia sobytia L
      t_ojz(L) = exprnd(1/lambda); % vremya, cherez kotoroe poyavitsya novaya zayavka, posle sobytia L
      j(L) = J += 1; % nomer zayavki j(L) = J

      J; % nomer zayavki J
      t_z(J) = t_sob(L); % moment poyavlenia zayavki J
      t_obsl(J) = exprnd(1/mu); % vremya obslujivania zayavki J
      t_kob(J) = t_z(J) + t_obsl(J); % moment okonchania obslujivania zayavki J
      % poisk svobodnogo pribora:
      for k(J) = 1:k0+1
        if apps(k(J)) == free
          apps(k(J)) = busy; %  -  zanyali pribor; zapomnili k(J)
          break
        endif
      endfor
      %  -  takim obrazom v k(J) zapishetsya nomer pribora, obslujivayuschego zayavku J)

      if ~isempty(doobsl)
        doobsl(2, :) -= t_ojz(L-1); %  -  vremena doobslujivania tekuschih zayavok umen'shayutsya na t_ojz!
      endif
      doobsl(:, end+1) = [J, t_obsl(J)]';
      [t_min(L), j_min] = min(doobsl(2, :)); %minimal'noe doobslujivanie
      freed = doobsl(1, j_min); % zapomnili nomer zayavki, obslujivanie kotoroy zakonchitsya ran'she ostalnyh

    case 2 % tip 2:

      t_sob(L) = t_sob(L-1) + t_min(L-1);
      t_ojz(L) = t_ojz(L-1) - t_min(L-1);
      j(L) = freed;
      apps( k(freed) ) = free; %  -  osvobodili sootvetstvuyuschy pribor
      if ~isempty(doobsl)
        doobsl(2, :) -= t_min(L-1); %  -  vremena doobslujivania vseh zayavok, ostavshihsya na obslujivanii, umen'shayutsya na t_min!
        [t_min(L), j_min] = min(doobsl(2, :));
        freed = doobsl(1, j_min);
        doobsl(:, find( doobsl(1, :) == freed)(end) ) = [];
      endif
  endswitch

  C(L) = sum(apps); % kolichestvo zayavok v obrabotke
  % minimal'noe ostavsheesya vremya obslujivania zayavok posle sobytia L:
  if C(L) == 0
    t_min(L) = -1;
  endif

  % obschee vremya prebyvania SMO v sostoyanii  i  s momenta t = 0 do t_sob(L):
  t(1:k0+1, L) = i = 0;
  while i < C(L-1)+1
    i += 1;
    t(i, L) = t(i, L-1);
  endwhile
  t(C(L-1)+1, L) = t_sob(L) - t_sob(L-1);

  delta = max(abs( t(:, L)' ./ t_sob(L) - r )); % otklonenie posle sobytia L
  delta_or_C = [delta, C(L)];


   %______________________________________used-in-used functions:
    function T = type(T_min, T_ojz)
      if T_ojz < T_min || T_min == -1
        T = 1;
      else
        T = 2;
      endif
    endfunction

  endfunction
end


function r = stationDist(po)
  r = [];
  err = 1;

  i = 0;
  while err >= 0.00001
    i+=1;
    r(i) = po^(i-1)/factorial(i-1)/exp(po);

    err = abs(r(i));
  endwhile

  r = [r, i]; % i = 0...k0, a u nas sdvinuto do 1...(k0+1) iz-za opredelenia massivov v oktave
endfunction


function print_it(arr)

  for i = 1:length(arr)
   printf('\n%i\n', round_it(arr(i), 6));
  endfor

  % okrugleniye do 5 znakov posle zapyatoy:
  function R = round_it(x, k)
    R = 0; R = round(x*10^k)/10^k;
  endfunction
end