function zust = zustandsbesetzer_Z1(MC,MK,PP,dt)

% Berechnet die Besetzung aller Zustände für äquidistante Zeitschritte
% C ist die C-Matrix noch ohne Konzentrationen und Hauptdiagonale
% MK - Konzentrationsprotokoll - ist die Matrix mit den Konzentrationen 
% PP ist das Pulsprotokoll
% dt ist der Zeitschritt

PoS = 0.985;

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);
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; 
end

% Besetzung der GGWB - Matrix für den ersten Puls
GGWB = zeros(ns,1);
CG = CC(:,:,1);
CG(end,:) = 1;
b = [zeros(ns-1,1);1];		% b-Vektor mit [0; 0; ...; 0; 1]
GGWB = CG\b;			% Zustandsvektor im GGW als Spalten-


% Zustandsbesetzung im 1. Puls im Gleichgewicht
nMP = PP(1,1)/dt + 1;
zust = ones(ns,nMP);
for i = 1:nMP
    zust(:,i) = GGWB;
end

% Zustandsbesetzung in den weiteren Pulsen
x0 = GGWB;          % Anfang 2. Puls entspricht Ende des 1. Pulses
for i = 2:np
    nMP = PP(1,i)/dt;
    CG = CC(:,:,i);
    t = (1:nMP)*dt;
    [ev,ew] = eig(CG);
    z = ev*(((ev\x0)*ones(1,nMP)).*exp(diag(ew)*t));
    zust = [zust,z];
    x0 = z(:,end);
end