close all
clear all
a = 1;
b = 1.5; % provare 1.5
% numero di simulazioni
J = 3000;
T = 50; % deve essere sufficientemente alto
f1 = 0;
tic
for j = 1:J
  t = 0;
% primo stato (0 persone in coda)
  s = 1;
  if (s == 1)
    b0 = 0;
  else
    b0 = b;
  end  
%  b0 = sign(s-1)*b;
  w = randexp(a+b0);
  t = t+w;
  while (t < T)
    s = s+2*randdisc([b0/(a+b0),a/(a+b0)])-3;
    if (s == 1)
      b0 = 0;
    else
      b0 = b;
    end  
%    b0 = sign(s-1)*b;
    w = randexp(a+b0);
    t = t+w;
  end  
  t = T;
  if (length (f1) < s)
    f1(s) = 1;
  else
    f1(s) = f1(s)+1;
  end
end
toc
smax1 = length(f1);
plot([1:smax1],f1(1:smax1)/J,'o')
% secondo modo
f2 = 0;
tic
for j = 1:J
  arrivi = zeros(ceil(2*a*T),1); % il tempo medio di interarrivo e` 1/a
  Narrivi = 1;
  lambdaa = randexp(a);
  arrivi(Narrivi) = lambdaa;
  while (arrivi(Narrivi) < T) 
    Narrivi = Narrivi+1;
    lambdaa = randexp(a);
    arrivi(Narrivi) = arrivi(Narrivi-1)+lambdaa;
  end  
  arrivi = arrivi(1:Narrivi);
  servizi = zeros(ceil(2*b*T),1); % il tempo medio di interservizio e` 1/b
  Nservizi = 1;
  lambdab = randexp(b);
  servizi(Nservizi) = lambdab;
  while (servizi(Nservizi) < T) 
    Nservizi = Nservizi+1;
    lambdab = randexp(b);
    servizi(Nservizi) = servizi(Nservizi-1)+lambdab;
  end  
  servizi = servizi(1:Nservizi);
% tempi di transizione
  tt = sortrows([[arrivi,ones(Narrivi,1)];[servizi,-ones(Nservizi,1)]]);
  s = 1;
  i = 0;
  t = tt(1,1);
  while (t < T)
    i = i+1;
    s = max(1,s+tt(i,2));
    t = tt(i+1,1);
  end
  t = T;
  if (length (f2) < s)
    f2(s) = 1;
  else
    f2(s) = f2(s)+1;
  end
end
toc
smax2 = length(f2);
smax = max(smax1,smax2);
hold on
rho = a/b;
plot([1:smax2],f2(1:smax2)/J,'*r',[1:smax],(1-rho)*rho.^[0:smax-1],'k')
legend('primo modo','secondo modo','invariante')