%Caso dei due servitori self

% Fissa i parametri di arrivo degli utenti e dei due server  e calcola tre serie aleatorie 
% con distribuzione esponenziale.

lambda_a=1/20; n0=100;
Ta=exprnd(1/lambda_a,n0,1);

lambda_s=1/11;
Ts1=exprnd(1/lambda_s,n0,1);
Ts2=exprnd(1/lambda_s,n0,1);


% Calcola i tempi di arrivo cumulati

ta(1)=Ta(1);

for k=2:n0
ta(k)=Ta(k)+ta(k-1);
end

%Calcola i tempi di servizio cumulati ai due server e il numero di utenti i coda presso i due server stessi

ts1(1)=ta(1)+Ts1(1);
ts2(1)=0; coda(1)=1; Nc1(1)=1; Nc2(1)=0; j=1; r=1; t=1; l=2;

for k=2:n0
ts1(k)=10000;
ts2(k)=10000;
end


   for k=2:n0

       while ta(k)>=ts1(j)
       j=j+1;
       end

       j=j-1;

       while ta(k)>=ts2(r)
       r=r+1;
       end

       r=r-1;

          if ts2(1)==0
          r=0;
          end

             coda(k)=k-(r+j);
 
m=1; n=1;

while ts1(m)<10000
       m=m+1;
end  

Nc1(k)=m-j-1;

while ts2(n)<10000
      n=n+1;
end  


Nc2(k)=n-r-1;


                if ts2(1)==0
                Nc2(k)=0;
                end

%A questo punto si invia il nuovo utente alla coda meno lunga nel case il numero di
%utenti in fila ai due server sia diverso. Nella situazione in cui ai due servitori
%ci siano invece lo stesso numero di utenti si attiva un processo aleatorio che invia
%l'ultimo arrivato ad uno dei due serventi con probabilità del 50%.

%Per ogni singolo caso viene calcolato il tempo cumulato di erogazione dei servizi .

if Nc1(k)<Nc2(k)
ts1(l)=max(ts1(l-1),ta(k))+Ts1(l);
l=l+1;
end

if Nc2(k)<Nc1(k)

    if t>1
    ts2(t)=max (ts2(t-1),ta(k))+Ts2(t);
    t=t+1;
    end

    if t==1
    ts2(1)=ta(k)+Ts2(1);
    t=t+1;
    end

end

if Nc1(k)==Nc2(k)

   if unidrnd(100)/100<=.5
   ts1(l)=max(ts1(l-1),ta(k))+Ts1(l);
   l=l+1;

   else

   if t>1
    ts2(t)=max (ts2(t-1),ta(k))+Ts2(t);
    t=t+1;
    end

    if t==1
    ts2(1)=ta(k)+Ts2(1);
    t=t+1;
    end

   end

end


j=1; r=1; n=1; m=1;

   end 


%Il codice seguente è applicabile, con i dovuti distinguo che saranno evidenziati nei commenti,
%sia al caso dei due serventi self sia a quello misto.

%Elimina dai tempi di servizio i valori 10000 assegnati precedentemente per ragioni
%di calcolo

for i=0:n0-1

   if ts1(n0-i)==10000
    ts1(n0-i)=[];
   end

   if ts2(n0-i)==10000
    ts2(n0-i)=[];
   end

end

%Nei casi 2 e 3 assegna i valori 1 ai tempi di arrivo e i valori -1 ai tempi di servizio.

TA=ta'; TS1=ts1'; TS2=ts2';

aux_vector=[TA;TS1;TS2];

ord_vector=sort(aux_vector);

for k=1:2*n0

ord_matrix(k,1)=ord_vector(k);

    for j=1:n0
    if ord_vector(k)==aux_vector(j)
    ord_matrix(k,2)=1;
    end
    end

     
    for j=n0+1:2*n0 
    if ord_vector(k)==aux_vector(j)
    ord_matrix(k,2)=-1;
    end
    end

end 

%Adesso effettua il calcolo degli utenti in coda istante per istante

for i=2:2*n0

   ord_matrix(i,2)=ord_matrix(i,2)+ord_matrix((i-1),2);

end

ord_matrix(:,1)=[];

%Comando facoltativo per plottare l'andamento della coda nel tempo

plot(ord_vector,ord_matrix)


%Calcola il numero di utenti serviti nelle due code nel tempo preso a riferimento. 
%Nel caso misto la cosa ha anche maggiore importanza dal momento che il ricavo 
%di una prestazione "servita" è maggiore.

ut_serv1=1; ut_serv2=1;

while TS1(ut_serv1)<180
ut_serv1=ut_serv1+1;
end

while TS2(ut_serv2)<180
ut_serv2=ut_serv2+1;
end


%Nel caso dei due serventi self il calcolo del ricavo si effettua semplicemente facendo:

ricavo=47*(ut_serv1+ut_serv2);