function strom = simulator_EN(N,f)

% Simuliert die Besetzung der Zustände in das Array zust für 1 Kanal

[MC,MK,uc,ub] = Modell_4Q;      % Modelldaten einlesen

[PP,dt] = f;     % hier möchte ich den Funktionsnamen übergeben lassen 
simulator_EN(5, @Puls_4Q);

AnzMP =    PP(1,:)/dt;          % Anzahl Messpunkte in Pulsen
AnzMPkum = cumsum(AnzMP);       % kumulative Summe

np = length(PP(1,:));           %  Anzahl Pulse
ns = length(MC(1,:));           % Anzahl Zustände

nMK = ones(ns,ns)-MK;           % Matrixwerte ohne Konzentrationen

% Besetzung des CC - Arrays
CC  = zeros(ns,ns,np);
ev3 = zeros(ns,ns,np);
ew3 = zeros(ns,np);

for k = 1:np
    
    konz = PP(3,k);             % Konzentration im Puls
    C = konz*MK.*MC+nMK.*MC;    % Multiplikation mit konz
    
    for i=1:ns                  % Hauptdiagonale von C berechnen
        C(i,i)=0;
        element=0;
        for j=1:ns
            element=element-C(j,i);
        end
        C(i,i)=element;
    end
    CC(:,:,k) = C;              % alle C-Matrizen
    [ev,ew]   = eig(C);         % Eigenvektoren und Eigenwerte
    ev3(:,:,k)= ev;            % np Matrizen mit Eigenvektoren
    ew3(:,k)  = diag(ew);      % np Spalten mit Eigenwerten
end


% Berechnet Gleichgewichtbesetzung für 1. Puls

GGWB1 = zeros(ns,1);
CG = CC(:,:,1);
CG(end,:) = 1;
b = [zeros(ns-1,1);1];          % b-Vektor mit [0; 0; ...; 0; 1]
GGWB1 = CG\b;                   % Zustandsvektor im GGW als Spaltenvektor

E = eye(ns);         % Einheitsmatrix, Spaltenvektoren sind Anfangszustände
D = triu(ones(ns));  % obere Dreiecksmatrix mit Einsen belegt
B = [];              % sicherheitshalber;

% Berechnet für jeden Puls und jeden möglichen Ausgangszustand (k)
% die Matrix B.pmat

for pu = 1:np
    EV = ev3(:,:,pu);    % Matrix der Eigenvektoren in Array
    EW = ew3(:,pu);      % Vektor der Eigenwerte in Matrix
        
    for k = 1:ns
        Anf_Besetzung = E(:,k);                      % nur Zustand k besetzt
        B(pu,k).pmat  = EV*(diag(EV\Anf_Besetzung)); % .pmat von typ und k
    end
end

% Berechnet für jeden trace (atr) und jeden Puls (pu) die summierten Wahr-
% scheinlichkeiten dafür, dass nach der Zeit dt der Zustand ns , 
% ns oder ns-1, ... , ns oder ns-1 oder ... oder 1 angenommen wird, wenn
% der Anfangszustand k war. Die summierten Wahrscheinlichkeiten stehen in 
% einer Matrix .sprobl, deren k-ter Spaltenvektor dem Ausgangszustand k 
% entspricht. 
% Abspeicherung der sprobls-Matrizen unter der Variablen tiS(atr,pu) 

for pu =1:np    
   EW = ew3(:,pu);
   for k = 1:ns                      % k ist ein bestimmter Anfangszustand
     probl = B(pu,k).pmat*exp(EW*dt); % Wahrscheinlichkeiten nach dt
     TIS(pu).sprobl(:,k) = D*probl;   % summierte Wahrscheinlichkeiten 
   end
end
  
% Die Vor-Besetzung ist ein Zustand, der entsprechend den Wahrscheinlich-
% keiten der GGWB eingenommen worden ist und aus dem heraus die Besetzung 
% des ersten Messpunktes im 1. Puls abgeleitet wird

zust_0 = zeros(ns,AnzMPkum(np));
zust_s = zust_0;                    % Summe der Einzelströme

for nr = 1:N
    zust = zust_0;                  % Nullmatrix für Zustände  
    sBesetzung = D*GGWB1;     % summierte Besetzungswahrscheinlichkeiten
    a = rand;                 % Zufallszahl entscheidet Vor - Besetzung   
    k = max(find(sBesetzung>a));    % zufälliger Zustand in Vor-Besetzung
    j = 0;
    for pu = 1:np         
       for i=1:AnzMP(pu)
                 j = j+1;
                 a=rand;
                 sBesetzung = TIS(pu).sprobl(:,k); 
                 % summierte Wahrscheinlichkeiten für Zustand in Punkt j
                 k = max(find(sBesetzung>a)); 
                 % k ist nun Nummer des in j besetzten Zustands
                 zust(:,j) = zust(:,j) + E(:,k);
                 % Spaltenvektor E(:,k) hat eine 1 an der k-ten Stelle
        end
    end
    zust_s = zust_s + zust;
end 
 
zust = zust_s / N;

strom = uc*zust;
% fluor = ub*zust;
t = dt*(1:length(strom));

plot(t,strom);