% Skript unwucht_5.m, in dem die Schwingungen durch Unwucht
% eines Feder-Masse-System mit Tilger untersucht werden
% Arbeitet mit Modell unwucht5.slx
clear;
s = tf('s');
% ------- Parameter des Systems
m = 10;   m0 = 1;
k = 10;   c = 5;    cT = 0.00001;
e = 0.2;
alpha = 0.1;      % Winkelbeschleunigung
% ------- Eigene Frequenz ohne Daempfung
omega_0 = sqrt(k/(m+m0));
f_0 = omega_0/(2*pi);
% ------- Frequenz der Anregung
%omegar = 0.5*omega_0;
%omegar = 2*omega_0;
omegar = omega_0;
fr = omegar/(2*pi);
% ------- Parameter des Tilgers (angepasst an omegar)
kT = 5,          
mT = kT/(omegar)^2,        % Angepasster Tilger an omegar
%mT = 0.9*kT/(omegar)^2,   % Nicht ganz angepasster Tilger
% ------- Frequenzbereich der Anregung
fmin = 0;       
fmax = fr;            % Stationäre Frequenz nach dem Anlauf
%fmax = 2.5*fr;        % Stationäre Frequenz nach dem Anlauf
% ------- 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;

[b1T,a1T] = tfdata(HT(1));
b1T = b1T{:};       % Koeffizienten der Uebertragungsfunktion 
a1T = a1T{:};       % von Fe bis x

[b2T,a2T] = tfdata(HT(2));
b2T = b2T{:};       % Koeffizienten der Uebertragungsfunktion 
a2T = a2T{:};       % von Fe bis xT

f = linspace(0,f_0*5,1000);
Hx = freqs(b1T, a1T, 2*pi*f);
figure(1);    clf;
subplot(211), plot(f, 20*log10(abs(Hx)));
title('Amplitudengang des Systems mit Tilger von Fe bis x');
xlabel('Hz');    ylabel('dB');     grid on;
axis tight;
subplot(212), plot(f, angle(Hx)*180/pi);
title('Amplitudengang des Systems mit Tilger von Fe bis x');
xlabel('Hz');    ylabel('Grad');     grid on;
axis tight;

% ------- Aufruf der Simulation
dt = 0.2;     Tfinal = 1000;
t = 0:dt:Tfinal;
T1 = 500;
alpha = 2*pi*(fmax-fmin)/T1,
f = (fmax - (fmax-fmin)*(t-T1)./(Tfinal-T1)).*(t>T1) +...
    fmax*(t<=T1);
my_options = simset('MaxStep', dt);
sim('unwucht5', [t], my_options);

t = x.time;
x  = x.signals.values;
xT = xT.signals.values;

figure(2);    clf;
subplot(311), plot(t, f);
title(['Frequenz der Anregung; Eigenfrequenz =',...
    num2str(f_0),' Hz']);
xlabel('Zeit in s');    grid on;    ylabel('Hz');
La = axis;     hold on;
plot(La(1:2), [f_0, f_0],'--');
hold off;

subplot(312), plot(t, x);
title('Schwingung der Lage der Hauptmasse');
xlabel('Zeit in s');    grid on;

subplot(313), plot(t, xT);
title('Schwingung der Lage der Tilgungsmasse');
xlabel('Zeit in s');    grid on;


















