% Aufgabe 29-12: Walze mit Feder gefesselt auf schiefer Ebene
%
% Die Function arbeitet mit den Standard-Einstellungen
%
%   Masse der Walze  m = 41 kg,
%   Integrationsintervall  t = 0 ... 20 s.
%
% Beide Werte koennen beim Aufruf (aus dem Command Window) veraendert
% werden, z. B.:
%
%   Aufg29_12Ani (30)    arbeitet mit  m = 30 kg  und  t = 0 ... 20 s.
%   Aufg29_12Ani (50,30) arbeitet mit  m = 50 kg  und  t = 0 ... 30 s.

function Aufg29_12Ani (mym , mytend) 

global m c alpha b R g 
c = 12000 ; alpha = pi/6 ; b = 0.8 ; R = 0.15 ; g = 9.81 ;

if nargin > 0  m    = mym    ;  else  m    = 43 ;  end
if nargin > 1  tend = mytend ;  else  tend = 20 ;  end

tanf = 0 ; xanf = 0 ; vanf =  0 ;                  % Anfangswerte

clf ;
options = odeset ('MaxStep'   , 0.01     , ...
                  'OutputFcn' , @odeplot , 'OutputSel' , [1]) ;
[t x]  = ode45 (@BewDglFederWalze , [tanf tend] , [xanf vanf] , options) ;

xi = x(:,1) ;
[t xi]
tt   = tanf ;
minx = -2*R ;
maxx = max(xi)*cos(alpha)+ 2*R  
miny = b/2 - max(xi)*sin(alpha)  
if (miny > 0)
    maxx = b/tan(alpha) 
    miny = -b/4 
end
maxy = b+1.5*R ;

clf ;
plot ([0 b/8 -b/8 0] ,[0 -b/6 -b/6 0], 'r-') ;    % Lager fuer die Feder
hold on
xe = -R*sin(alpha)+max(xi)*cos(alpha) ;
ye = b-R*cos(alpha)-max(xi)*sin(alpha) ;
if (ye > 0)
    ye = 0 ;
    xe = (b-R*cos(alpha))/tan(alpha)-R*sin(alpha) ;
end
plot ([-R*sin(alpha) xe] ,[b-R*cos(alpha) ye], 'r-') ;  % Schiefe Ebene
hold on

ncirc = 24 ;                      % Anzahl der Polygonpunkte fuer Kreis
[xs ys nc] = SpringCoords ([0 0] ,[xi(1) b] , 4 , 2*b) ; % Federkoordinaten
xs = zeros(nc+ncirc+1) ;
[xs(nc+1:nc+ncirc+1) ys(nc+1:nc+ncirc+1) nc2] = ...
             CylCoords (xi(1) , b , R , ncirc , 0) ;    % Walzenkoordinaten  
p  = plot (xs , ys , 'EraseMode' , 'xor') ;
h  = findobj ;
axis equal ; axis ([minx maxx miny maxy]) ;             % Zeichenfeldlimit
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 (t) ;

for ii=1:Zeitschritte    % Schleife über alle berechneten Zeitschritte
   if (t(ii) >= tt)      % Ausgabe jeweils nur nach deltatt 
       [xs ys nc] = SpringCoords ([0 0] ,...
                 [xi(ii)*cos(alpha) b-xi(ii)*sin(alpha)] ,  4 , 2*b) ;
       [xs(nc+1:nc+ncirc+1) ys(nc+1:nc+ncirc+1) nc2] = ...
             CylCoords (xi(ii)*cos(alpha) , b-xi(ii)*sin(alpha) , ...
                        R , ncirc , xi(ii)/R) ;
       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 xvp = BewDglFederWalze (t , xv)

global m c alpha b R g

x = xv(1) ;                % Input-Werte kommen in
v = xv(2) ;                % einem Vektor an.

Lc = sqrt(b^2+x^2-2*b*x*sin(alpha)) ;  % Aktuelle Federlaenge
cb = (x-b*sin(alpha))/Lc ;             % cos (beta)
Fc = c*(Lc-b) ;                        % Federkraft

xp = v ;                               % Die beiden Differenzialgleichungen
vp = (m*g*sin(alpha)-Fc*cb)/(1.5*m) ;  % werden ausgewertet, ...

xvp = [xp ; vp] ;  % ... die Ergebnisse werden in einem Vektor abgeliefert.