%% Der umfallende Schornstein

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016
% Version 2.5  2011


close all, clear, clc
%% Materialwerte
rho = 7.86;                     % kg/dm^3 Eisen
l = 10;                         % m Laenge
r = 10/2;                       % cm Radius
%
m = rho * pi * (r/10).^2 *l*10; % Masse
J = 1/3 * m * l.^2;             % Traegheitsmoment
g = 9.81;                       % m/s^2

%% Ende Erdoberflaeche
tmax = 5.;%30.;%8.;
y0=[1*pi/180;0];%[30*pi/180;0];%[10*eps;0];%[pi/180/60;0]; ohne kleine Auslenkung bleibt der Schornstein stehen
figure
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',0.1, ... 
    'Events',@(t,y) kippStabevent(t,y,l,m,J));
[t,y] = ode23(@(t,y) kippStabDGL(t,y,l,m,J),[0:0.005:tmax],y0,options);
% aequidistante Zeitschritte wegen Visualisierung

%% Visualisierung 
he = l * cos(y(1,1)) - 1/2 * g * t.^2;
h2e = l * cos(y(:,1));
plot(t,he,t,h2e),shg
xlabel('Zeit'),ylabel('Hoehe')
legend('freies Teilchen','Stab','location','NorthWest')
title('Vergleich Freier Fall-Stab')
%
figure,plot(t,h2e), grid on
xlabel('Zeit'),ylabel('Hoehe')
ve = y(:,2)/J*l;                  % Geschwindigkeit Stab
vef = sqrt(2*l*cos(y(1,1))*g);    % max. Geschwindigkeit freier Fall
vee = sqrt(3*l*cos(y(1,1))*g);    % max. Geschwindigkeit Stab
vteste = ve - vef;
figure, plot(vteste,h2e),grid on
xlabel('Geschwindigkeit - Maxmimum'),ylabel('Hoehe')
title('Differenzgeschwindigkeit')

%% Veranschaulichung der Beschleunigung
figure
nlauf = 1:2:length(t);
n=max(nlauf)-100;%650;
    % Plotten des Stabs
    xp = l * linspace(0,sin(y(n,1)),25);
    yp = l * linspace(0,cos(y(n,1)),25);
    plot(xp,yp),axis equal, xlim([0,10]),ylim([0,10]),shg, hold on
    % Plotten der Beschleunigung
    lp = sqrt(xp.^2+yp.^2);      % Position auf dem Stab
    ap = lp/l * 3/2 * g * sin(y(n,1));  % Beschleunigung ddphi
    up = ap * cos(y(n,1));              % x Komponente
    vp = -ap * sin(y(n,1));             % y Komponente
    xp(end+1) = l;                      % dient der Skalierung
    yp(end+1) = 0;
    up(end+1) = 0;
    vp(end+1) = -sqrt(3*l*g);
    quiver(xp,yp,up,vp),shg
    axis equal
    %
    ug = g * cos(y(n,1))*ones(size(xp));
    vg = -g * sin(y(n,1))*ones(size(xp));
    ug(end) = 0;
    vg(end) = vp(end);
    quiver(xp,yp,ug,vg)
    %
    ugg = up-ug;
    vgg = vp-vg;
    ugg(end) = 0;
    vgg(end) = vp(end);
    quiver(xp,yp,ugg,vgg)
    title('Veranschaulichung der Kraefte')