BeginPackage["queue3`"]
queue3::usage="This is a collection of functions used in this book."
gammadist::usage="Computes the value of the gamma distribution using
the notation in this book."
hyperg::usage="The pmf of the hypergeometric distribution."
hypergdist::usage="The distribution function of the hypergeometric
distribution."
mmc::usage="Computes the performance statistics for the M/M/c queueing system."
mmck::usage="Computes the performance statistics for the M/M/c/K queueing
system."
bca::usage="Erlang's B (loss) formula."
cca::usage="Erlang's C (delay) formula."
cent::usage="Computes the performance statistics for the central server
model with fixed MPL N."
online::usage="Computes the performance statistics for terminal system
with FESC to replace the central server model of the computer system.
subcent can be used to calculate the rates needed as input."
subcent::usage="Computes the throughput for a central server model with
fixed MPL."
Begin["queue3`private`"]
gammadist/: gammadist[alpha_, beta_, t_] :=
Gamma[beta, 0, alpha*t]/Gamma[beta]

hyperg[n_, nsucc_, ntot_, k_] :=
Binomial[nsucc, k]*Binomial[ntot - nsucc, n - k]/Binomial[ntot, n]

hypergdist[n_, nsucc_, ntot_, x_] :=
Sum[hyperg[n, nsucc, ntot, k], {k, 0, Min[x, nsucc, ntot]}]

bca[c_, a_]:=
Block[{t=1,n},
Do[t= a t /(n+ a t), {n, 1, c}];
Return[t]
]
cca[c_,a_]:=
Block[ {n, b, cca},
rho =a/c ;
b=a/(1+a);
For[n=2, n<=c, n++, b= a*b/(n+a*b)] ;
cca=b/(rho*b+1-rho);
Return[{cca}] ;
]

mmc[ c_, x_, s_ ] :=
(* x is lambda, s is mean service time *)
Block[ {n, a, w, rho, b, cca, QueueTime },
a = x s ;
rho =a/c ;
b=a/(1+a);
For[n=2, n<=c, n++, b= a*b/(n+a*b)] ;
cca=b/(rho*b+1-rho);
QueueTime = cca*s/(c*(1-rho)) ;
w=QueueTime+s;
(* QueueTime is mean time in the queue *)
(* w is mean time in the system *)
Return[{QueueTime, w, rho, cca}] ;
]

mmck[ c_, k_, x_, s_ ] :=
Block[ {cf, kf, n, a, p, p0, pk, q, Xa, QueueTime },
cf = Floor[c] ;
kf = Floor[k] ;
a = x s ;
p = 1 ;
p0 = 1 ;
q = 0 ;
For[n=1, n<=cf, n++, p = p a/n ; p0 = p0 + p ] ;
For[n=c+1, n<=kf, n++, p = p a/cf ; p0 = p0 + p ; q = q + (n-cf) * p ] ;
pk = p / p0 ;
Xa = x * (1-pk) ;
QueueTime = (q / p0 ) / Xa ;
w=QueueTime + s;
(* pk is the probability a customer is lost *)
(* Xa is true throughput *)
(* QueueTime is mean time in queue *)
(* w is mean time in system *)
Return[{pk, Xa, QueueTime, w}] ;
]

cent[k_,N_, D_]:=
(* central server model *)
(* k is number of service centers *)
(* N is MPL, D is service demand vector *)
Block[{L, w, wn, n, lambdan, rho},
L=Table[0, {k}];
For[n=1, n<=N, n++, w=D*(L+1); wn=Apply[Plus,w]; lambdan=n/wn;
L=lambdan w; rho=lambdan D];
(* lambdan is mean throughput *)
(* wn is mean time in system *)
(* L is vector of number at servers *)
(* rho is vector of utilizations *)
Return[{lambdan, wn, L, rho}];
]

online/: online[srate_, N_, T_]:=
Block[{n, w,s},
m=Length[srate];
x=Table[Last[srate], {N-m}];
nsrate=Join[srate, x];
q=Join[{1}, Table[0, {N-1}]];
s=0;
q0=1;
For[n=1, n<=N, n++,
w=0;
For[j=1, j<=n, j++,
w=w+(j/nsrate[[j]])*If[j>1, q[[j-1]], q0];
lambda=n/(T+w)];
s=0;
For[j=n, j>=1, j--,
q[[j]]=(lambda/nsrate[[j]])*If[j>1, q[[j-1]],q0];
s=s+q[[j]]];
q0=1-s
];
qplus=Join[{q0},q];
(* lambda is mean throughput *)
(* w is mean response time *)
(* qplus is vector of conditional probabilities *)
Return[{lambda,w, qplus}];
]

subcent/: subcent[k_,N_, D_]:=
Block[{L, w, wn, n, lambdan, rho},
L=Table[0, {k}];
For[n=1, n<=N, n++, w=D*(L+1); wn=Apply[Plus,w]; lambdan=n/wn;
L=lambdan w; rho=lambdan D];
Return[{lambdan}];
]

End[] (* end 'private' context *)
EndPackage[]