% Skript foerder_haft_gleit_2.m in dem die Förderung mit 
% Haft- und Gleitreibung simuliert wird. Arbeitet mit 
% Modell foerder_haft_gleit1.mdl
clear;
% -------- Parameter des Systems
g = 9.89;
% ------- Anregungsbewegung
ampl = 0.02;     % Amplitude der Anregungsbewegung
fanr = 4;      % Frequenz der Anregungsbewegung
mu = 0.5;        % Reibungskoeffizient
% ------- Simulation 
Tfinal = 5;
dt = 0.0001;
t = 0:dt:Tfinal;    nt = length(t);
% Initialisierungen
z = ampl*cos(2*pi*fanr*t);
vz = -ampl*2*pi*fanr*sin(2*pi*fanr*t);
az = -ampl*(2*pi*fanr)^2*cos(2*pi*fanr*t);
beta = pi/4;
x = z*cos(beta);      y = z*sin(beta);
vx = vz*cos(beta);    vy = vz*sin(beta);
ax = az*cos(beta);    ay = az*sin(beta);
% ------- Signale für das Modell
ay_m = [t', ay'];       vx_m = [t', vx'];
% ------- Aufruf der Simulation
Tstart = 0;      Zeitschritt = dt;
simOut = sim('foerder_haft_gleit2','StartTime','Tstart',...
    'StopTime','Tfinal','Solver','ode2','FixedStep',...
    'Zeitschritt','SaveOutput','on','OutputSaveName','yout_new',...
    'SaveTime','on','TimeSaveName','tout_new');
y = simOut.get('y');    % 'y' Name der Variablen aus To Workspace
% Ausgaänge auf Outport:
yout_new = simOut.get('yout_new'); % Ergibt eine Matrix mit drei Spalten
% Man muss im Menü 'Model Configuration Parameters' die Option 'Limit
% data points to last' deaktivieren
tout_new = simOut.get('tout_new'); % Simulationszeit 
%######################
t = y.time;
vxx = y.signals.values(:,1);
xx = y.signals.values(:,2);
Fr = y.signals.values(:,3);

figure(1);   clf;
nperioden = 5;
Tnperioden = 5/fanr;     nT = fix(Tnperioden/dt);
nd = nt - nT:nt;         % Darstellungen von 5 Perioden im
                         % stationären Zustand
% nd = 1:nt;             % Darstellung mit Einschwingen
axen(1) = subplot(411), plot(t(nd),(g+ay(nd)));
title(['Summe der Beschleunigungen (g+ay)']);
xlabel('Zeit in s');   grid on;   axis tight;
axen(2) = subplot(412), plot(t(nd), vxx(nd), t(nd), vx(nd),...
    t(nd),10*x(nd));
title(['Geschwindigkeiten vxx und vx' ,...
    ' bzw. Lage des Förderorgans 10*x']);
xlabel('Zeit in s');   grid on;   axis tight;
axen(3) = subplot(413), plot(t(nd), Fr(nd),...
    t(nd),sign(vxx(nd)-vx(nd)'));
title('Reibungskraft / m und Signum(vxx - vx)');
xlabel('Zeit in s');   grid on;   axis tight;
axen(4) = subplot(414), plot(t(nd), xx(nd));
title('Bewegung des Förderguts');
xlabel('Zeit in s');   grid on;   axis tight;
linkaxes([axen(1),axen(2),axen(3),axen(4)],'x');


