% Feder-Masse-System mit Reibung.

function Aufg28_5Ani (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.8 ;  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 ...
             & xve(telen,1) > a              % 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  myq=my ;     % Ereignis "Ab sofort Reibung"
      elseif ereig(telen) == 3  myq=0  ;     % Ereignis "Ab sofort reibungsfrei" 
      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

x = xvt(1) ;
v = xvt(2) ;
value(1)      =  v  ; % Ereignis 1: "Geschwindigkeit Null"
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

% ===================================================================
% 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] ;

% ===================================================================
function [x , y , nc] = SpringCoords (xyA , xyB , n , sl)

% ... berechnet Polygon (Koordinaten werden in x und y 
% abgeliefert, nc ist die Anzahl der Polygonpunkte) eines
% Federsymbols zwischen den Punkten xyA und xyB. n ist
% die Anzahl der gewuenschten "Zacken" (n=4 bedeutet 4
% "Zacken" auf jeder Seite, also viermal "von der Mitte
% nach aussen, von dort zur anderen Seite nach außen und
% von dort zurueck zur Mitte). sl ist die komplette
% "Drahtlaenge der Feder".

dist = sqrt ((xyB(1)-xyA(1))^2+(xyB(2)-xyA(2))^2) ;

nc = 2*n+4 ;
dx = (xyB(1)-xyA(1))/nc ;
dy = (xyB(2)-xyA(2))/nc ;
ds = dist/nc ;
slq = sl-3*ds ;
ql  = slq/(2*n) ;
b   = .5*sqrt(ql^2-ds^2) ;

x(1)    = xyA(1) ;
x(nc)   = xyB(1) ;
y(1)    = xyA(2) ;
y(nc)   = xyB(2) ;
x(2)    = xyA(1)+dx*1.5 ;
x(nc-1) = xyB(1)-dx*1.5 ;
y(2)    = xyA(2)+dy*1.5 ;
y(nc-1) = xyB(2)-dy*1.5 ;

ca = (xyB(1)-xyA(1))/dist ;
sa = (xyB(2)-xyA(2))/dist ;
xx = x(2) + dx ;
yy = y(2) + dy ;

for i=3:2:nc-3
    x(i) = xx - b*sa ;
    y(i) = yy + b*ca ;
    xx   = xx + dx ;
    yy   = yy + dy ;
    x(i+1) = xx + b*sa ;
    y(i+1) = yy - b*ca ;
    xx   = xx + dx ;
    yy   = yy + dy ;
end

% ===================================================================
function [x , y , nc] = MassCoords (xyC , b , h)

% ... berechnet Polygon (Koordinaten werden in x und y 
% abgeliefert, nc ist die Anzahl der Polygonpunkte) eines
% Rechtecks mit der Breite b und der Hoehe h, dessen 
% Mittelpunkt durch die Koordinaten xyC bestimmt ist.
% Das Polygon startet im Mittelpunkt.

b2   = b/2 ;
h2   = h/2 ;
x(1) = xyC(1) ;
y(1) = xyC(2) ;
x(2) = xyC(1) ;
y(2) = xyC(2)+h2 ;
x(3) = xyC(1)+b2 ;
y(3) = xyC(2)+h2 ;
x(4) = xyC(1)+b2 ;
y(4) = xyC(2)-h2 ;
x(5) = xyC(1)-b2 ;
y(5) = xyC(2)-h2 ;
x(6) = xyC(1)-b2 ;
y(6) = xyC(2)+h2 ;
x(7) = xyC(1) ;
y(7) = xyC(2)+h2 ;
x(8) = xyC(1) ;
y(8) = xyC(2)-h2 ;
nc   = 8 ;