% Feder-Masse-System mit Reibung.

function FederMasseReibungAni2 (myin , my0in)

% Parameter (global definiert, damit Wertaenderung nur an
% einer Stelle erforderlich):
global m c my a g r myq
m = 20 ; c = 300 ; a = 0.5 ; g = 9.81 ; 

if nargin > 0  my  = myin  ;  else  my  = 0.4 ;  end
if nargin > 1  my0 = my0in ;  else  my0 = my  ;  end
tend = 10 ;

x0    = [2*a ; 0] ;          % Anfangswerte
tspan = [0 tend]    ;        % Zeitintervall

if x0(1) > 0    r = -1 ; else   r = 1 ; end
if x0(1) > a  myq = my ; else myq = 0 ; end

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

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

while cont == 1
  [t xv te xve ereig]  = ode45 (@FederMasseReibung_Dgl , 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) ; 
      tspan = [te(telen) tend] ;             % ... weiter mit "Restintervall ...
      x0    = [xve(telen,1) xve(telen,2)] ;  % ... und neuen Anfangsbedingungen
      if ereig(telen) == 1                   % Ereignis "Geschwindigkeit Null"
          if abs(c*xve(telen,1)) < my0*m*g   % Federkraft kleiner Haftkraft, ...
              cont = 0 ;                     % ... also Abbruch
              tout  = [tout  ; tend]  ;
              xvout = [xvout ; [xv(lent,1) 0]] ;
          elseif xve(telen,1) < 0  r=1 ;     % Es geht nach rechts weiter ...
          else  r=-1 ;                       % ... bzw. nach links
          end
      elseif ereig(telen) == 2               % Ereignis "Punkt x=a erreicht"
          if xve(telen,2) > 0  myq=my ;  else  myq=0 ;  end
      end
  end
end

% "Sortieren der Ergebnisse":
xi = [1 0] * xvout' ;  % xvout' ist die transponierte Matrix der Ergebnisse,

tt = tout(1) ;
maxi = max(xi) ; if (maxi < 0) maxi =  a/4 ; end
mini = min(xi) ; if (mini > 0) mini = -a/4 ; end
clf ;
plot ([mini-a mini-a max(xi)+0.5*a] ,[0.3*a -0.3*a -0.3*a], 'r-') ;
hold on
plot ([a max(xi)+0.5*a] , [-0.3*a -0.3*a], 'g-') ;
hold on
[xs ys nc] = SpringCoords ([mini-a 0] ,[xi(1) 0] ,  6 , 5*a) ;
[xs(nc+1:nc+8) ys(nc+1:nc+8) nc2] = MassCoords ([xi(1) 0] , 0.8*a , 0.5*a) ;
p  = plot (xs , ys , 'EraseMode' , 'xor') ;
h  = findobj ;
axis ([1.2*mini-a 1.2*maxi+a -a/4 a/4]) ; axis equal ;
zeit = title (sprintf('t = %8.2f s' , tt)) ;
set (gca , 'userdata' , zeit)

deltatt  = 1/25     ;    % 25 Bilder pro Sekunde
zeitlupe = 1        ;    % Verzögerungsfaktor
cstart   = cputime  ;
Zeitschritte = length (tout) ;

for ii=1:Zeitschritte    % Schleife über alle berechneten Zeitschritte
   if (tout(ii) >= tt)      % Ausgabe jeweils nur nach deltatt 
       [xs ys nc] = SpringCoords ([mini-a 0] ,[xi(ii) 0] ,  6 , 5*a) ;
       [xs(nc+1:nc+8) ys(nc+1:nc+8) nc2] = MassCoords ([xi(ii) 0] , 0.8*a , 0.5*a) ;
       set(p , 'XData' , xs , 'YData' , ys)
       while (cputime-cstart < tt*zeitlupe)  % Warteschleife
       end 
       zeit = get (gca , 'userdata');
       set(zeit , 'string' , sprintf('t = %8.2f s' , tt))
       drawnow
       tt = tt + deltatt ;
   end
end

% ===================================================================
function [value,isterminal,direction] = Ereignis (t,xvt)
global m c my a g r myq

value(1) = xvt(2) ; % Ereignis 1: "Geschwindigkeit Null"
isterminal(1) = 1 ; % Abbruch, ...
direction(1)  = 0 ; % ... in jedem Fall 

value(2) = xvt(1)-a ; % Ereignis 2: "Punkt x=a erreicht"
isterminal(2) = 1 ;   % Abbruch ...
direction(2)  = 0 ;   % in jedem Fall

% ===================================================================
% Funktion, die die "rechte Seite" des Dgl.-Systems definiert: 
function f = FederMasseReibung_Dgl (t , xvt)

% Input: Aktueller Zeitpunkt t und die Werte fuer Weg xi und 
% Geschwindigkeit vi zu diesem Zeitpunkt im Vektor  xvt = [xi ; vi]

% Output: xPunkt und vPunkt im Vektor  f = [xPunkt ; vPunkt] 

global m c my a g r myq

xi = xvt(1) ;
vi = xvt(2) ;

xPunkt = vi ;
vPunkt = - c*xi/m - r*myq*g ;

f = [xPunkt ; vPunkt] ;
