function p=controlQ(seed,mu)
% estimates proportion of time that an MG1 queue >20 in system using MM1 queue 
% for control variate sampling. arrival rate=1, service time is gamma(.9,1)
% mu=rate for eponential service times.
IA=exprnd(1,1,1000);
U=rand(1,1000);
S1=gaminv(U,1.1,1);
S2=-(1/mu)*log(1-U);
[ANS,ATS,arrivet,departt,ET1,NS1]=MG1Q(IA,S1);
T1=Time(20,ET1,NS1);
[ANS,ATS,arrivet,departt,ET2,NS2]=MG1Q(IA,S2);
T2=Time(20,ET2,NS2);
% the expected number in the MM1 system is rho/(1-rho)=1/(mu-1)
% and the probability of more than 20 is sum{(1-rho)rho^n;n=21.....}=(1-rho)^21
p=T1-T2+(1-1/mu)^21;   % this is a control variate estimator