% Skript aktiv_tilger_12.m, in dem eine aktive Tilgung
% untersucht wird. Arbeitet mit Modell aktiv_tilger12.mdl
clear;
s = tf('s');
% -------- Parameter des Systems
m1 = 200;      k1 = 100;     c1 = 20;
f1 = sqrt(k1/m1)/(2*pi),    % Eigenfrequenz ohne Dämpfung
m2 = 1;       k2 = 10;      c2 = 10;
f2 = sqrt(k2/m2)/(2*pi),    % Eigenfrequenz ohne Dämpfung
fsig = 0.1;     % Frequenz der Anregung Fe1
Fe1 = 0.5;     % Amplitude der Anregung
% -------- Regler mit Rückführung von der Beschleunigung a1 
P = 0;       I = 2000;     D = 0; % Mit Regler (aktive Tilgung)
%P = 0;       I = 0;     D = 0;     % Ohne Regler 
% -------- Vereinfachte Übertragungsfunktionen
H1 = P + I/s + D*s/(s/100+1);     
H2 = (c2*s+k2)/(m2*s^2+c2*s+k2);  
H3 = s^2/(m1*s^2+c1*s+k1);

% -------- Übertragungsfunktion der geöffneten Schleife
Hoff = H1*H2*H3;
Hoff = minreal(Hoff);
[boff, aoff] = tfdata(Hoff);
boff = boff{:},    aoff = aoff{:},
roots(aoff),
% Frequenzgang
f = logspace(-2,2);     w = 2*pi*f;
Hoffs = freqs(boff, aoff, w);

figure(1);   clf;
subplot(211), semilogx(f, 20*log10(abs(Hoffs)));
title('Amplitudengang der geöffneten Schleife');
xlabel('Hz');     ylabel('dB');     grid on;

subplot(212), semilogx(f, angle(Hoffs)*180/pi);
title('Phasengang der geöffneten Schleife');
xlabel('Hz');     ylabel('Grad');     grid on;
% -------- Übertragungsfunktion der geschlossenen Schleife
Hges = H1*H2*H3/(1+H1*H2*H3);
Hges = minreal(Hges);
[bges, ages] = tfdata(Hges);
bges = bges{:},    ages = ages{:},
roots(ages),
% Frequenzgang
f = logspace(-2,2);     w = 2*pi*f;
Hges = freqs(bges, ages, w);

figure(2);   clf;
subplot(211), semilogx(f, 20*log10(abs(Hges)));
title('Amplitudengang der geschlossenen Schleife');
xlabel('Hz');     ylabel('dB');     grid on;

subplot(212), semilogx(f, angle(Hges)*180/pi);
title('Phasengang der geschlossenen Schleife');
xlabel('Hz');     ylabel('Grad');     grid on;

% -------- Aufruf der Simulation
Tfinal = 1000;     dt = 0.01;
sim('aktiv_tilger12', [0:dt:Tfinal]);
t = y.time;
y1 = y.signals.values(:,1);
y2 = y.signals.values(:,2);
Fe1_t = f.signals.values(:,1);  % Anregungskraft
Fe2_t = f.signals.values(:,2);

% -------- Für die Darstellung von 10 Perioden 
%          im stationären Zustand
nt = length(t);         nd = 1:nt; 
Tdarst = 10/fsig;       ndarst = fix(Tdarst/dt);
if ndarst <= nt
     nd = nt - ndarst : nt;
else nd = 1:nt;
end;
Q = std(y1(nd)),      % Gutefunktion
figure(3);    clf;
subplot(311), plot(t(nd),y1(nd));
  title(['Koordinate y1(t) für fsig = ',num2str(fsig), ' Hz'])
  xlabel('Zeit in s');      grid on;
subplot(312), plot(t(nd),y2(nd));
  title(['Koordinate y2(t)'])
  xlabel('Zeit in s');      grid on;
subplot(313), plot(t(nd),y1(nd),t(nd),y2(nd));
  title(['Koordinate y1(t) und y2(t)'])
  xlabel('Zeit in s');      grid on;
 
figure(4);    clf;
plot(t(nd),Fe1_t(nd),t(nd),Fe2_t(nd));
  title(['Anregungs- und Rückführungskraft (Fe1, Fe2) für fsig = ',...
      num2str(fsig), ' Hz']);
  xlabel('Zeit in s');      grid on;
  legend('Fe1', 'Fe2')

% -------- Übertragungsfunktion von der Störung bis zur Lage y1
Hstoe = H3/((1+H1*H2*H3)*s^2);
Hstoe = minreal(Hstoe);
[bstoe, astoe] = tfdata(Hstoe);
bstoe = bstoe{:},    astoe = astoe{:},
roots(astoe),
% Frequenzgang
f = logspace(-2,2);     w = 2*pi*f;
Hstoes = freqs(bstoe, astoe, w);

figure(5);   clf;
subplot(211), semilogx(f, 20*log10(abs(Hstoes)));
title(['Amplitudengang von der Störung Fe1' ...
    ' bis zur Koordinate y1']);
xlabel('Hz');     ylabel('dB');     grid on;

subplot(212), semilogx(f, angle(Hstoes)*180/pi);
title('Phasengang von der Störung bis zur Beschleunigung');
xlabel('Hz');     ylabel('Grad');     grid on;

