# Simulating a (G,G,1) queue
n <- 2000 # number of customers
T <- rgamma(n,shape=1,scale=1) #interarrival times
s <- 0.8 #scale of service times (mean if shape =1)
S <- rgamma(n,shape=1,scale=s) #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
l <- 10 # start counting at l-th event
w <- E[(l+1):k] - E[l:(k-1)] # inter-event times between l-th and k-th event
N.r <- N[(l+1):k] # customers during these intervals
M <- max(N.r) # maximal number of customers
p <- rep(0,M+1)
for (j in (0:M)) p[j+1] <- sum(w*(N.r==j))/(E[k]-E[l])
plot(0:M,p,type="h",xlab="k",ylab="P[ k customers]",
     main="Distribution of number of customers")
#comparison with analytical result for shape = 1
points(0.1+0:M,(1-s)*s^(0:M),type="h",col=2)