# Simulating a (G,G,1) queue
n <- 1000 # number of customers
T <- rgamma(n,shape=1,scale=1) #interarrival times
S <- rgamma(n,shape=1,scale=0.7) #service times

A <- cumsum(T) #arrival times
W <- rep(0,n) # workload at times A[i]-
for (i in (2:n)) W[i] <- max(0,W[i-1]+S[i-1]-T[i])
D <- A+S+W #departure times
N <- c(rep(1,n),rep(-1,n)) #changes in the number of customers at A and D
E <- c(A,D) #event times (arrivals or departures)
ord <- order(E)
E <- E[ord] #ordered event times
N <- c(0,cumsum(N[ord])) # number of customers at event times

#Plots
#number of customers
f.N <- stepfun(E,N)
plot(f.N,do.points=FALSE,xlab="time",ylab="number of customers",
     main="Queueing system")
#workload
plot(A,W+S,xlab="time",ylab="workload")
points(A,W)
abline(h=0)
segments(A,W,A,W+S)
segments(A,W+S,D,rep(0,n))

#Statistics
#number of customers
k <- (1:(2*n))[E==A[n]] # last arrival = k-th event
w <- E[1:k] - c(0,E[1:(k-1)]) # inter-event times until last arrival
N <- N[1:k]
M <- max(N) # maximal number of customers
p <- rep(0,M+1)
for (j in (0:M)) p[j+1] <- sum(w*(N==j))/A[n]
plot(0:M,p,type="h",main="Distribution of number of customers")