function [S,Q,a,tt] = codagg1(T)
%
% [S,Q,a] = codagg1(T)
%
% interarrivo come gamma(alphasigma)/lambdasigma
alphasigma = 1;
lambdasigma = 1.4;
a = 1/(alphasigma/lambdasigma);
% servizio come gamma(alphatau)/lambdatau
alphatau = 1;
lambdatau = 2;
b = 1/(alphatau/lambdatau);
rho = a/b; % funziona bene per rho vicino a 1
n = 1;
t = 0;
sigma(n) = randg(alphasigma)/lambdasigma;
A(n) = sigma(n);
while (A(n) < T)
  n = n+1;
  sigma(n) = randg(alphasigma)/lambdasigma;
  A(n) = A(n-1)+sigma(n);
end
A = A(1:n-1);
nclienti = n-1;
tau = randg(alphatau*ones(1,nclienti))/lambdatau;
W(1) = 0;
for n = 1:nclienti-1
  W(n+1) = max(W(n)+tau(n)-sigma(n+1),0);
end
S = W+tau;
arrivi = [A(1:nclienti)',ones(nclienti,1)];
uscite = [A(1:nclienti)'+S',-ones(nclienti,1)];
% tempi in cui succede qualcosa
tt = sortrows([arrivi;uscite]);
countermax = length(tt);
counter = 1;
Q(counter) = tt(counter,2); % Q(1)=1
counter = counter+1;
t = tt(counter,1);
% i tempi di uscita possono superare T:  bisogna fermarsi prima
while ((t < T) && (counter < countermax))
  Q(counter) = Q(counter-1)+tt(counter,2);
  counter = counter+1;
  t = tt(counter,1);
end
t = T;
tt = tt(:,1)';
tt = tt(1:length(Q));
