%
%      MAUT.M
%
%  Verhalten einer Mautstation
%
% Jan Lunze
% 3.4.2006
% für 2. Auflage: 6.7.2012
% für 3. Auflage: 26.1.2017
% für 4. Auflage: 12.1.2025
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  lambda = Parameter des Ankunftsprozesses
%  mu = Parameter des Bedienprozesses
clear
close all
echo on
%
%   Bedien- und Ankunftsprozesse
%
echo off
%  Ausnutzung der Eigenschaft, dass die Zeitpunkte, an denen einer
%  Poissonprozess seinen Zustand ändert, über das betrachtete Zeitintervall
%  gleichverteilt sind
lambda=6/60;    % 6 Fahrzeuge pro Minute
Nmax=100;
Tend=Nmax/lambda;
load('RandMaut');
rng(RandMaut);
TA=Tend*rand(1, Nmax);
TA=sort(TA);
TA=[0 TA];
Z=[0:Nmax];
for k=1:Nmax-1
    TDelta(k)=TA(k+1)-TA(k);
end
mu=8/60;         % 8 Fahrzeuge pro Minute
Tend=Nmax/mu;
TC=Tend*rand(1,Nmax);
TC=sort(TC);
TC=[0 TC];
Z=[0:Nmax];
for k=1:Nmax-1
    TB(k)=TC(k+1)-TC(k);
end
N=0;
T=0;
TW=0;
for k=1:Nmax-1
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
end
TD=TA(1:Nmax-1)+TW(1:Nmax-1)+TB(1:Nmax-1);
figure(1)
subplot(7,1,[6,7])
stem(TW, 'LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0]);
%    mittlere Wartezeit
TWquer=sum(TW)/length(TW);
Tmax=max(TD);
plot([0 Nmax+1],[TWquer, TWquer],':');
latexylabel('$$T_{\rm W}$$');
rotateY;
yticks([0 25 50]);
latexxlabel('$$k$$');
axis([0 100 -0.5 max(TW)+0.5]);
hold off
%
%
subplot(7,1,[1,2])
queueplot(TA, TD);
xleer;
%latexxlabel('$$t$$ in s');
xticksempty([0 200 400 600 800 1000]);
axis([0 1000 -0.3 0.3]);
latextitle('Wartesystem mit stochastischen Zeiten ($$\lambda=$$0,100, $$\mu=$$0,133)');
%
%  Bestimmung der Anzahl N von Kunden im Wartesystem
%
kA=0;
kD=0;
kAmax=length(TA);
kDmax=length(TD);
kmax=kAmax+kDmax;
N(1)=0;
T(1)=0;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=TD(kD+1);
        kD=kD+1;
    end
end       % while-Schleife verlassen, wenn entweder alle TA-Werte 
          % oder alle TD-Werte berücksichtigt sind
while kD<kDmax    % Auswertung der restlichen TD-Werte
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TA(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N, 'b', 'LineWidth',1.5);
hold on
%   mittlere Anzahl von Kunden im System
Nquer=tmean(N, T);
plot([0 Tmax],[Nquer, Nquer],'b--');
plot([0 Tmax],[0 0]);
axis([0 1000 -0.5 max(N)+0.3]);
latexylabel('$$N$$');
rotateY;
yticks([0 6]);
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off;
%   Überprüfung des Gesetzes von Little
NW=max(N-1,0);
NWquer=tmean(NW, T)
TWquer

%%
epsfigc('Maut');