% Skript unwucht_2.m, in dem die Schwingungen durch Unwucht
% eines Feder-Masse-System mit Tilger untersucht werden
clear;
s = tf('s');
% ------- Parameter des Systems
m = 10;   m0 = 1;
k = 10;   c = 5;    cT = 0.001;
e = 0.2;
% ------- Eigene Frequenz ohne Daempfung
omega_0 = sqrt(k/(m+m0));
f_0 = omega_0/(2*pi);
% ------- Frequenz der Anregung
%omega = 0.5*omega_0;
%omega = 5*omega_0;
omega = omega_0;
f = omega/(2*pi);
% ------- Parameter des Tilgers (angepasst an omega)
kT = 5,          
mT = kT/(omega)^2,      % Angepasster Tilger
mT = 0.9*kT/(omega)^2,   % Nicht ganz angepasster Tilger
% ------- Amplitude Hauptmasse bei omega ohne Tilger
H1 = 1/((m+m0)*s^2+c*s+k);
[b1,a1] = tfdata(H1);
b1 = b1{:};       % Koeffizienten der Uebertragungsfunktion
a1 = a1{:};
zaehler1 = polyval(b1, j*omega);
nenner1  = polyval(a1, j*omega);
ampl_x = m0*e*omega^2*abs(zaehler1/nenner1), % Ampl. ohne Tilger

% ------- Uebertragungsfunktionen mit Tilger
% Matrizen des Systems
Ai = [(m+m0)*s^2+(c+cT)*s+(k+kT), -cT*s-kT; -cT*s-kT, mT*s^2+cT*s+kT];
Bi = [1, 0]';
HT = inv(Ai)*Bi;
% ------- Amplitude Hauptmasse x bei omega mit Tilger
[b1T,a1T] = tfdata(HT(1));
b1T = b1T{:};       % Koeffizienten der Uebertragungsfunktion 
a1T = a1T{:};       % von Fe zu x
zaehler1T = polyval(b1T, j*omega);
nenner1T  = polyval(a1T, j*omega);
ampl_xT = m0*e*omega^2*abs(zaehler1T/nenner1T), % Ampl. mit Tilger
% ------- Amplitude Tilger xT bei omega 
[b2T,a2T] = tfdata(HT(2));
b2T = b2T{:};       % Koeffizienten der Uebertragungsfunktion 
a2T = a2T{:};       % von Fe zu xT
zaehler2T = polyval(b2T, j*omega);
nenner2T  = polyval(a2T, j*omega);
ampl_xmT = m0*e*omega^2*abs(zaehler2T/nenner2T), % Ampl. des Tilgers

gewinn = ampl_x/ampl_xT,    % Gewinn
% -------- Amplitudenfunktion von x ohne Tilger abhängig von omega
f_t = linspace(f/10, f*10, 10000);
omega_t = 2*pi*f_t;
zaehler1 = polyval(b1, j*omega_t);
nenner1  = polyval(a1, j*omega_t);
ampl_x = m0*e*omega_t.^2.*abs(zaehler1./nenner1);
figure(1);  
subplot(211), plot(f_t, ampl_x);
title('Amplitude x ohne Tilger');
xlabel('Hz');     grid on;
hold on;     La = axis;
plot([f, f], [La(3), La(4)],'r');   hold off;
subplot(212), plot(f_t, ampl_x);
title('Amplitude x ohne Tilger (Ausschnitt)');
xlabel('Hz');     grid on;
hold on;     La = axis;
plot([f, f], [La(3), La(4)],'r');   hold off;
axis([f*0.7, f*1.3, La(3:4)]);
% -------- Amplitudenfunktion von x mit Tilger abhängig von omega
zaehler1T = polyval(b1T, j*omega_t);
nenner1T  = polyval(a1T, j*omega_t);
ampl_xT = m0*e*omega_t.^2.*abs(zaehler1T./nenner1T);
figure(2);  
subplot(211), plot(f_t, ampl_xT);
title('Amplitude x mit Tilger');
xlabel('Hz');     grid on;
hold on;     La = axis;
plot([f, f], [La(3), La(4)],'r');   hold off;
subplot(212), plot(f_t, ampl_xT);
title('Amplitude x mit Tilger (Ausschnitt)');
xlabel('Hz');     grid on;
hold on;     La = axis;
plot([f, f], [La(3), La(4)],'r');   hold off;
axis([f*0.7, f*1.3, La(3:4)]);
% -------- Amplitudenfunktion von xT des Tilgers abhängig von omega
zaehler2T = polyval(b2T, j*omega_t);
nenner2T  = polyval(a2T, j*omega_t);
ampl_xmT = m0*e*omega_t.^2.*abs(zaehler2T./nenner2T);
figure(3);  
subplot(211), plot(f_t, ampl_xmT);
title('Amplitude xT des Tilgers');
xlabel('Hz');     grid on;
hold on;     La = axis;
plot([f, f], [La(3), La(4)],'r');   hold off;
subplot(212), plot(f_t, ampl_xmT);
title('Amplitude xT des Tilgers (Ausschnitt)');
xlabel('Hz');     grid on;
hold on;     La = axis;
plot([f, f], [La(3), La(4)],'r');   hold off;
axis([f*0.7, f*1.3, La(3:4)]);










