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

function Aufg29_12 (mym , myxanf) 

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

if nargin > 0  m    = mym    ;  else  m    = 50 ;  end
if nargin > 1  xanf = myxanf ;  else  xanf =  0 ;  end

tanf = 0 ; tend = 10 ; vanf =  0 ;      % Integrationsintervall, Anfangswert
     
% Runge-Kutta-Algorithmus -------------------------------------------------- 
nRK = 1000                              % Anzahl der Runge-Kutta-Schritte
[t xv] = rk4 (@BewDglFederWalze , [tanf tend] , [xanf vanf] , nRK) ;
% Runge-Kutta-Algorithmus --------------------------------------------------

% MATLAB-ode45-Solver mit MaxStep-Option -----------------------------------
maxstep  = 0.01  
options  = odeset ('MaxStep' , maxstep) ;
[to xvo] = ode45 (@BewDglFederWalze , [tanf tend] , [xanf vanf] , options) ;
% MATLAB-ode45-Solver mit MaxStep-Option -----------------------------------

x  = xv (:,1) ;                           % Sortieren der Ergebnisse:
v  = xv (:,2) ;                           % Spalte 1: x, Spalte 2: v
xo = xvo(:,1) ;                           % Sortieren der Ergebnisse:
vo = xvo(:,2) ;                           % Spalte 1: xo, Spalte 2: vo

xEndRK    = x(end) 
vEndRK    = v(end) 
xEndode45 = xo(end) 
vEndode45 = vo(end) 

Lc        = sqrt(b^2+x.^2-2*b*x*sin(alpha)) ;
TgesRK    = 0.75*m*v.^2 - m*g*x*sin(alpha) + 0.5*c*(Lc-b).^2 ;
Lcode45   = sqrt(b^2+xo.^2-2*b*xo*sin(alpha)) ;
Tgesode45 = 0.75*m*vo.^2 - m*g*xo*sin(alpha) + 0.5*c*(Lcode45-b).^2 ;

subplot(3,1,1) , plot(t , x , to , xo , '-r') ; grid on ,
title('Weg-Zeit-Gesetz') , ylabel('x') , xlabel('t') ; 
subplot(3,1,2) , plot(t , v , to , vo , '-r') ; grid on ,
title('Geschwindigkeits-Zeit-Gesetz') , ylabel('v') , xlabel('t') ;
subplot(3,1,3) ,
plot(t , TgesRK , to , Tgesode45 , '-r') ;
title('Gesamtenergie') , ylabel('Tges') , xlabel('t') ;
legend ('Runge-Kutta 4. Ordnung', 'MATLAB-ode45' , 3) ;

% =========================================================================
function xvp = BewDglFederWalze (t , xv)

global m c alpha b 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)) ;
cb = (x-b*sin(alpha))/Lc ;
Fc = c*(Lc-b) ;

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.