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

clear;
% -------- Parameter des Systems
Fex_ampl = 1;
m1 = 100;        k1 = 5000;     rg = 0.5;
w01 = sqrt(k1/m1);
m2 = 5;          mu = m2/m1;
k2 = mu*k1/((1+mu)^2);
c2 = 2*m2*sqrt(k2/m2)*sqrt(3*mu/(8*(1+mu)));  % optimale Dämpfung
v0 = 0.0;        y0 = 0.0;
v02 = 0;         y02 = 0
% -------- Aufruf der Simulation
Tfinal = 500;
dt = 0.001;      ts = [0,Tfinal];
my_options = simset('solver','ode4','FixedStep',dt);
%my_options = simset('solver','ode45','MaxStep',dt);

sim('tilger_gleitreib_21',ts,my_options);

y1  = y.signals.values(:,1);
y2 = y.signals.values(:,2);
t = y.time;

figure(1);   clf; %############
a1 = subplot(211), plot(t,y1);
title('Koordinate der Lage der Masse 1');
xlabel('Zeit in s');     grid on;
a2 = subplot(212), plot(t,y2);
title('Koordinate der Lage der Tilgermasse 2');
xlabel('Zeit in s');     grid on;
linkaxes([a1,a2],'x');
Q = sum(y1.^2),     % Gütekriterium
% -------- Spektrale Leistungsdichte
nfft = 256;
[Pss1,f] = pwelch(y1(1:50:end),hamming(nfft),50,nfft,1/(50*dt));
[Pss2,f] = pwelch(y2(1:50:end),hamming(nfft),50,nfft,1/(50*dt));
figure(2);    clf;
plot(f, 10*log10([Pss1, Pss2]));
title(' Spektrale Leistungsdichte in dB/Hz');
xlabel('Hz');    ylabel('10*log10(Pss)');  grid on;
legend('Masse 1','Masse 2');
