%
%  macm 202 -- 04 feb 03 -- djm
%
%  w6gap1.m:  simulate street crossing process
%
%     - use direct exponential distribution
%

clear

%  simulation of Nruns

Nruns = 10000;
mu = 1.0;  K = 1;  walk = K/mu;

%  bins for histogram & pseudorandom seed

bins = (0:1/8/mu:5*walk);
rand('state',7654321);

wait_times = zeros(Nruns,1);
car_data = zeros(Nruns,1);

%  simulation loop

for k=1:Nruns;
    wait = 0;
    cars = 0;
    gap = -log(rand(1,1))/mu;
    while (gap<walk)
        wait = wait + gap;
        cars = cars + 1;
        gap = - log(rand(1,1))/mu;
    end
    wait_times(k,1) = wait;
    car_data(k,1) = cars;
end    

%  plot histogram of queue_size probabilities

figure(1); clf
count = hist(wait_times,bins);
bar(bins,count/Nruns,'hist')
hold on

title('\bf simulation of street crossing')
xlabel('\bf waiting time')
ylabel(['\bf probability for wait time = ' num2str(walk*mu) '\times 1/\mu'])
axis([0 max(bins) 0 1.1*max(count/Nruns)])

out = [mean(wait_times) var(wait_times)]

%  car data

figure(2); clf
cbins = 0:24;
count = hist(car_data,cbins);
bar(cbins,count/Nruns,'hist','b')
hold on

title('\bf simulation of street crossing')
xlabel('\bf # of cars')
ylabel(['\bf probability of cars for walk time = ' num2str(walk*mu) '\times 1/\mu'])
axis([0 max(cbins) 0 1.1*max(count/Nruns)])

plot(cbins,exp(-K)*(1-exp(-K)).^cbins,'r*')

out = [mean(car_data) var(car_data)]

theory = [exp(K)-1 ; (exp(K)-1)^2 + (exp(K)-1)]

%  car data

figure(3); clf
plot(cbins,log(count/Nruns+10^(-16)),'bo')
hold on

title('\bf simulation of street crossing')
xlabel('\bf # of cars')
ylabel(['\bf log(probability of cars)'])
axis([0 max(cbins) -10 0])

plot(cbins,log(exp(-K)*(1-exp(-K)).^cbins),'r*')
hold on
slope = log(1-exp(-K));
plot([0 15],[0 15]*slope-0.5)
text(5,-2,['\bf slope = ' num2str(slope)])