function LaufkatzeEvKontr ()

% Parameter (global definiert, damit Wertaenderung nur an einer Stelle erforderlich):
global mK mL JL ct lS Dt Ft a g
mK = 100 ; mL = 500 ; JL = 400 ; c = 200000 ; lS = 4 ; Dt = 1 ; F0 = 2000 ; a = 5 ; g = 9.81 ;

x0    = [0 ; 0 ; 0 ; 0] ;  ct = 0 ; Ft = F0 ;       % Anfangswerte
tend  = 10 ;      
tspan = [0 tend] ;                                  % Zeitintervall

% Integration des Anfangswertproblems:
options = odeset ('Mass' , @Massenmatrix , 'MaxStep' , 0.01 , 'Events' , @Ereignis) ; 

tout = [] ; xvout = [] ; % Auf tout und xvout werden die Ergebnisse gesammelt
cont = 1 ;

while cont == 1 
  [t xv te xve ereig]  = ode45 (@RechteSeite , tspan , x0 , options) ; 
  tout  = [tout  ; t]  ; % Die gerade abgelieferten Ergebnisse werden an ...
  xvout = [xvout ; xv] ; % ... die bereits gespeicherten Ergebnisse angehaengt

  lent = length (t) ;
  if (t(lent) >= tend)
     cont = 0 ;          % Ende erreicht, Abbruch der while-Schleife
  else                   % ... war ein "Ereignis" der Grund fuer den Abbruch
     telen = length (te) ; 
     disp(['Ereignis ' , num2str(ereig(telen)) , '  nach  t = ' , ... 
            num2str(te(telen)) , ' s  bei  x = ' , num2str(xve(telen)) , ' m']) ;

     if ereig(telen) == 1                   % Ereignis "Eine Sekunde nach dem Start" 
        Ft = 0 ;                            % Kraft abschalten
     elseif ereig(telen) == 2 ct = c ;      % ... Feder zuschalten
     elseif ereig(telen) == 3 ct = 0 ;      % ... Feder abschalten 
     end

     tspan = [te(telen) tend] ;             % ... weiter mit "Restintervall ...
     x0   = [xve(telen,1) xve(telen,2) xve(telen,3) xve(telen,4)] ; 
  end                                       % ... und neuen Anfangsbedingungen
end

% Grafische Ausgabe von x und phi, zunaechst "Sortieren der Ergebnisse":
x     = xvout(:,1) ;  % xvout ist die Matrix der Ergebnisse,
phi   = xvout(:,2) ;  % x, phi ... 
v     = xvout(:,3) ;  % v und omega ... 
omega = xvout(:,4) ;  % ... sind Vektoren 

Tges  = 0.5*mK*v.^2 + 0.5*mL*(v.^2+lS^2*omega.^2+2*lS*v.*omega.*cos(phi)) ...
       +0.5*JL*omega.^2-mL*g*lS*cos(phi) ;           % Kontrollfunktion "Gesamtenergie"   

subplot(3,1,1) ; plot (tout , x)    , grid on , title  ('Bewegung der Laufkatze: x(t)') , axis([0 10 -25 10])
subplot(3,1,2) ; plot (tout , phi)  , grid on , title  ('Pendeln der Last: \phi(t)')
subplot(3,1,3) ; plot (tout , Tges) , grid on , title  ('Gesamtenergie: T_{ges}(t)')
disp(['Ergebnisse nach ' , num2str(length(tout)-1) , ' Zeitschritten:'])  ; % "End"-Kontrolle
disp(['x    (t_end = ' , num2str(tout(end)) , ') = ' , num2str(x  (end))])  ; 
disp(['phi  (t_end = ' , num2str(tout(end)) , ') = ' , num2str(phi(end))])  ;
disp(['T_ges(t_end = ' , num2str(tout(end)) , ') = ' , num2str(Tges(end))])  ;
disp(['(zum Vergleich: T_ges(t = ' , num2str(0) , ') = ' , num2str(Tges(1)) , ')'])  ;

% ===================================================================
% Function, die die Ereignisse definiert: 
function [value,isterminal,direction] = Ereignis (t,xv) 
global mK mL JL ct lS Dt Ft a g

x = xv(1) ;            % Aktuelle Position der Laufkatze

value(1)      = t-Dt ; % Ereignis 1: "Eine Sekunde nach dem Start"
isterminal(1) =   1  ; % Abbruch, ...
direction(1)  =   0  ; % ... in jedem Fall 

value(2)      = x-a ;  % Ereignis 2: "Punkt x=a erreicht"
isterminal(2) =   1 ;  % Abbruch, ...
direction(2)  =   1 ;  % ... wenn Bewegung nach rechts

value(3)      = x-a ;  % Ereignis 3: "Punkt x=a erreicht"
isterminal(3) =   1 ;  % Abbruch, ...
direction(3)  =  -1 ;  % ... wenn Bewegung nach links

% ============================================================
% Function, die die "Massenmatrix" des Dgl.-Systems definiert: 
function M = Massenmatrix (t , xv)

global mK mL JL ct lS Dt Ft a g

phi = xv(2) ;

M   = [1  0       0                0        ; 
       0  1       0                0        ; 
       0  0     mK+mL        mL*lS*cos(phi) ; 
       0  0  mL*lS*cos(phi)    mL*lS^2+JL   ] ;

% =================================================================== 
% Function, die die "rechte Seite" des Dgl.-Systems definiert: 
function r = RechteSeite (t , xv)

global mK mL JL ct lS Dt Ft a g


x = xv(1) ; phi = xv(2) ; v = xv(3) ; omega = xv(4) ;

b1 = mL*lS*omega^2*sin(phi)-ct*(x-a)+Ft ;
b2 = -mL*g*lS*sin(phi) ;

r = [v ; omega ; b1 ; b2] ;