% Skript aktiv_tilger_3.m, in dem eine aktive Tilgung
% untersucht wird. Arbeitet mit Modell aktiv_tilger3.mdl
clear;
s = tf('s');
% -------- Parameter des Systems
m1 = 200;      k1 = 100;     c1 = 20;
f1 = sqrt(k1/m1)/(2*pi),
m2 = 1;       k2 = 10;      c2 = 10;
f2 = sqrt(k2/m2)/(2*pi),
fsig = 0.1;
% -------- Analyse des Systems
Ai = [m1*s^2+(c1+c2)*s + (k1+k2), -(c2*s+k2);
    -(c2*s+k2), m2*s^2+c2*s+k2];
Bi = [1,0;0,1];

H = inv(Ai)*Bi;
H = minreal(H);
[b11, a11] = tfdata(H(1,1));   [b12, a12] = tfdata(H(1,2));
[b21, a21] = tfdata(H(2,1));   [b22, a22] = tfdata(H(2,2));
b11 = b11{:};    a11 = a11{:};
b12 = b12{:};    a12 = a12{:};
b21 = b21{:};    a21 = a21{:};
b22 = b22{:};    a22 = a22{:};
% ------- PID-Regler
P = 0;    I = 0;     D = 2000;
%P = 0;    I = 0;     D = 0;
Hsys = H(1,1)*feedback(1, (P+I/s+D*s/(s/100+1))*H(1,2)); 
% Übertragungsfunktion von der Störung Fe1 bis zur Lage y1
Hsys = minreal(Hsys);           
[bsys, asys] = tfdata(Hsys);
bsys = bsys{:},    asys = asys{:},
roots(asys),
% -------- Aufruf der Simulation
Tfinal = 1000,         dt = 0.01;
sim('aktiv_tilger3', [0:dt:Tfinal]);
t = y.time;
y11 = y.signals.values(:,1);
y12 = y.signals.values(:,2);
y1 = y.signals.values(:,3);
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)),    % Gütefaktor
figure(1);    clf;
subplot(311), plot(t(nd),y11(nd));
title(['Koordinate y11  ( fsig = ',num2str(fsig),' Hz )']);
xlabel('Zeit in s');    grid on;   axis tight;
subplot(312), plot(t(nd),y11(nd),t(nd),y12(nd));
title('Koordinate y11 und y12');
xlabel('Zeit in s');    grid on;   axis tight;
subplot(313), plot(t(nd),y1(nd));
title('Koordinate y1 = y11+y12');
xlabel('Zeit in s');    grid on;   axis tight;
figure(2);    clf;
f = logspace(-3,2,500);
Hs = freqs(bsys, asys, 2*pi*f);
subplot(211), semilogx(f, 20*log10(abs(Hs)));
title(['Amplitudengang von der Störung Fe1 bis zur Koordinate y1'...
    ' fsig = ', num2str(fsig), ' Hz und f2 = ', num2str(f2) ' Hz'])
xlabel('Hz');    ylabel('dB');    grid on;   %axis tight;
subplot(212), semilogx(f, angle(Hs)*180/pi);
title(['Phasengang von der Störung Fe1 bis zur Koordinate y1'])
xlabel('Hz');    ylabel('Grad');  grid on;   %axis tight;
% -------- Vereinfachte Übertragungsfunktion
H1 = D*s/(s/100+1);
H2 = (c2*s + k2)/(m2*s^2 + c2*s + k2);
H3 = 1/(m1*s^2+c1*s+k1);

Hsys1 = feedback(H3, H1*H2);      % Übertragungsfunktion von der
Hsys1 = minreal(Hsys1);           % Störung Fe1 bis zur Lage y1
[bsys1, asys1] = tfdata(Hsys1);
bsys1 = bsys1{:};    asys1 = asys1{:};

figure(3);     clf;
Hs1 = freqs(bsys1, asys1, 2*pi*f);
subplot(211), semilogx(f, 20*log10(abs(Hs1)));
hold on;
subplot(211), semilogx(f, 20*log10(abs(Hs)),'r');
title(['Amplitudengänge von der Störung Fe1 bis zur Koordinate y1 ('...
  ' Eigenfreq. f1 = ', num2str(f1), ' Hz und f2 = ', num2str(f2) ' Hz)'])
xlabel('Hz');    ylabel('dB');    grid on;   %axis tight;
hold off;
subplot(212), semilogx(f, angle(Hs1)*180/pi);
hold on;
subplot(212), semilogx(f, angle(Hs)*180/pi,'r');
title(['Phasengänge von der Störung Fe1 bis zur Koordinate y1'])
xlabel('Hz');    ylabel('Grad');  grid on;   %axis tight;
hold off;



