clear;

max_time = 50;
lambda = 3;
mu = 5;

interarrivals = exprnd(1/lambda, 1, floor(max_time*lambda*1.1));
arrivals(1)=interarrivals(1);
for j=2:length(interarrivals)
    arrivals(j) = arrivals(j-1) + interarrivals(j);
end
interdepartures = exprnd(1/mu, 1, floor(max_time*mu*1.1));
departures(1)=interdepartures(1);
for j=2:length(interdepartures)
    departures(j) = departures(j-1) + interdepartures(j);
end

arrivals(2,:) = 1;
departures(2,:) = -1;
c = [arrivals departures];
[a,idx] = sort(c, 2);
events = zeros(size(c));
for j=1:length(c)
    events(:,j)=c(:,idx(1,j));
end

maxidxs = find(events(1,:)>max_time);
maxidx = maxidxs(1)-1;
events(3,1)=0;
for j=2:maxidx
    events(3,j)=max(0, events(3,j-1)+events(2,j));
end

subplot(2,1,1);
plot(events(1,1:maxidx),events(3,1:maxidx));
subplot(2,1,2);
[h,x]=hist(events(3,1:maxidx));
plot(x,h/maxidx);

avglen = mean(events(3,1:maxidx))