%%%%% Coda ad un servente.
%%% Arrivi secondo Poisson di intertempo medio 2 gg
%%% quindi parametro lambda=1/2.
%%% Tempo di lavorazione esponenziale di media 1.5 gg
%%% quindi parametro mu=1/1.5.

m.a <-  2
n0 <- 100
T.a <- rweibull(n0,1,m.a)
hist (T.a, 20)

%%% hist (T.a, 20) per sincerarsi della correttezza.

m.s <- 1.5
T.s <- rweibull(n0,1,m.s)
hist (T.s, 20)


%%% Il tempo medio di servizio è inferiore a quello di interarrivo, quindi
%%% il sistema dovrebbe stabilizzarsi.

%%% Istanti di arrivo:

t.a <- matrix(nrow=n0)
for (k in 1:n0) {
t.a[k] <- sum(T.a[1:k])
} 
plot(t.a)


%%% Idea della risoluzione:
t.s[1] <- t.a[1]+T.s[1] è l'istante del primo servizio;
t.s[2] <- max(t.s[1], t.a[2])  + T.s[2]  è l'istante del secondo servizio;
t.s[3] <- max(t.s[2], t.a[3])  + T.s[3]  è l'istante del terzo servizio;
%%% e così via.

t.s <- matrix(nrow=n0)
t.s[1] <- t.a[1]+T.s[1]
for (k in 2:n0) {
t.s[k] <- max(t.s[k-1], t.a[k])  + T.s[k]
} 
plot(t.s)

%%% Numero di utenti in coda agli istanti di arrivo (immediatamente dopo ogni arrivo):

Na <- matrix(nrow=n0)

%%% Idea: quando arriva il k-esimo cliente, quanti dei precedenti (k-1)
sono stati serviti? tanti quanti i numeri t.s[j] < t.a[k].

Na[1] <- 1
for (k in 2:n0) {
Na[k] <- k - sum((sign(t.a[k]-t.s)+1)/2)
}
plot(Na)

%%% Elenco più sistematico: ordiniamo gli eventi temporali specificando 
per ciascuno di essi se si tratta di un arrivo o di una partenza:

S <- matrix(nrow=2*n0,ncol=2)

%%% Porremo S[k,1] = istante dell'evento k-esimo, 
%%% S[k,2] = 1 se è un arrivo, -1 se è un servizio.

S[,1] <- sort(c(t.a,t.s))


for (k in 1:(2*n0)) {
S[k,2]<-sum( (-1)*(sign(-abs(t.s-S[k,1])+0.001)+1)/2+1*(sign(-abs(t.a-S[k,1])+0.001)+1)/2)
}

plot(S[,1],S[,2])


%%% Con questi possiamo tracciare il grafico del numero di utenti nel sistema
ad ogni istante:

N <- matrix(nrow=2*n0)
for (k in 1:(2*n0)) {
N[k] <- sum(S[1:k,2])
}
plot(S[,1],N)

%%% Calcoliamo alcune frequenze sperimentalmente: 

pi.emp <- matrix(nrow=1000)
DS <- S[2:(2*n0),1]-S[1:(2*n0-1),1]
for (k in 1:1000) {
pi.emp[k] <- sum(DS*(sign(-abs(k-N[1:(2*n0-1)])+0.001)+1)/2)/S[(2*n0),1]
}


%%%%%
%%%%% Tempo di lavorazione Weibull di shape=5, scale=m.s

m.a <-  5
n0 <- 100
T.a <- rweibull(n0,1,m.a)
m.s <- 1.5
T.s <- rweibull(n0,5,m.s)
hist (T.s, 20)

t.a <- matrix(nrow=n0)
for (k in 1:n0) {
t.a[k] <- sum(T.a[1:k])
} 
t.s <- matrix(nrow=n0)
t.s[1] <- t.a[1]+T.s[1]
for (k in 2:n0) {
t.s[k] <- max(t.s[k-1], t.a[k])  + T.s[k]
} 
S <- matrix(nrow=2*n0,ncol=2)
S[,1] <- sort(c(t.a,t.s))
for (k in 1:(2*n0)) {
S[k,2]<-sum( (-1)*(sign(-abs(t.s-S[k,1])+0.001)+1)/2+1*(sign(-abs(t.a-S[k,1])+0.001)+1)/2)
}
N <- matrix(nrow=2*n0)
for (k in 1:(2*n0)) {
N[k] <- sum(S[1:k,2])
}
plot(S[,1],N)


%%%%%
%%%%% Tempo di lavorazione gaussiano di sd=0.01, media=m.s


m.a <-  5
n0 <- 100
T.a <- rweibull(n0,1,m.a)
m.s <- 1.5
T.s <- rnorm(n0,m.s,0.01)
hist (T.s, 20)

t.a <- matrix(nrow=n0)
for (k in 1:n0) {
t.a[k] <- sum(T.a[1:k])
} 
t.s <- matrix(nrow=n0)
t.s[1] <- t.a[1]+T.s[1]
for (k in 2:n0) {
t.s[k] <- max(t.s[k-1], t.a[k])  + T.s[k]
} 
S <- matrix(nrow=2*n0,ncol=2)
S[,1] <- sort(c(t.a,t.s))
for (k in 1:(2*n0)) {
S[k,2]<-sum( (-1)*(sign(-abs(t.s-S[k,1])+0.001)+1)/2+1*(sign(-abs(t.a-S[k,1])+0.001)+1)/2)
}
N <- matrix(nrow=2*n0)
for (k in 1:(2*n0)) {
N[k] <- sum(S[1:k,2])
}
plot(S[,1],N)


%%% si vede che il sistema è molto più "deterministico".