% Skript tilger_gleitreib3.m, in dem ein Tilger
% mit Gleitreibung untersucht wird.
% Arbeitet mit Modell tilger_gleitreib_3.mdl

clear;
% -------- Parameter des Systems
m = 10;     
k = 10;     
rg = 8,      m1 = 0.1;    k1 = 5,
%rg = 0;     k1 = 0,      % Ohne Isolierer
f0 = sqrt(k/m)/(2*pi),
f = f0;
ampl = 0.01;              % Für Sinusanregung
v0 = 0;      x0 = 0.002;
% -------- Aufruf er Simulation
Tfinal = 1000;       dt = 0.001;
my_options = simset('solver','ode4','MaxStep', dt);
sim('tilger_gleitreib_3',[0,Tfinal], my_options);
x = y(:,1);
x1 = y(:,2);
nx = length(x);
nd = 1:fix(nx/2);      % Bereich für die Darstellung          

figure(1);   clf;
subplot(211), plot(t(nd),x(nd));
title('Koordinate der Hauptmasse')
xlabel('Zeit in s');   grid on;
subplot(212), plot(t(nd),x1(nd));
title('Koordinate der Zwischenmasse')
xlabel('Zeit in s');   grid on;

% ------- Zielfunktion
Q = mean(detrend(x).^2),  % Gütekriterium
Q1 = std(x)^2,

% ------- Spektralanalyse
nf = fix(nx/2);    % Einschwingen entfernen
window = 16*4096;
h = spectrum.welch('Hamming',window, 50);

X = psd(h, x(nf:end),'Fs',1/dt);
X1 = psd(h, x1(nf:end), 'Fs',1/dt);
Xf = X.data;             Xf1 = X1.data;
f = X.frequencies;
nd = fix(window/200);     % Bereich für die Darstellung
figure(2),    clf;
subplot(211), plot(f(1:nd), 10*log10(Xf(1:nd)));
title('Spektrale Leistungsdichte der Hauptmasse');
xlabel('Hz');    ylabel('dB/Hz');     grid on;
La = axis;   axis([La(1:2), -150, La(4)]);
subplot(212), plot(f(1:nd), 20*log10(Xf1(1:nd)));
title('Spektrale Leistungsdichte der Zwischenmasse');
xlabel('Hz');    ylabel('dB/Hz');     grid on;
La = axis;   axis([La(1:2), La(3:4)]);
