% Skript vibration_control_1.m in dem eine Tilgung
% mit elektromagnetischem Dämpfer untersucht wird
clear;
s = tf('s');

% ------- Parameter des Systems
m1 = 100;    k1 = 5000;     c1 = 20;
m2 = 5;      mu = m2/m1;
k2 = mu*k1/((1+mu)^2),   % Nach Den Hartog
kf = 10;      kv = kf;
w01 = sqrt(k1/m1),      w02 = sqrt(k2/m2),
f01 = w01/(2*pi);
L = 0.5;    R = 15;    
% ------- Übertragungsfunktionen mit RL-Netzwerk
Ai = [m1*s^2+c1*s+(k1+k2),-k2,kf;
    -k2, m2*s^2+k2,-kf; 
    -kv*s, kv*s, L*s+R];
Bi = [1;0;0];
H = inv(Ai)*Bi;   H1 = minreal(H);
H11 = H1(1);      H12 = H1(2);     H13 = H1(3);
% ------- Koeffizienten der Übertragungsfunktionen
[b11, a11] = tfdata(H11);      [b12, a12] = tfdata(H12);
[b13, a13] = tfdata(H13);
b11 = b11{:};       a11 = a11{:};  % Y1/Fex  
b12 = b12{:};       a12 = a12{:};  % Y2/Fex
b13 = b13{:};       a13 = a13{:};  % I/Fex
% ------- Pole der Übertragungsfunktion Y1/Fex
roots(a11),
% -------- Klassische Tilgung zum Vergleich
% Parameter nach Den Hartog
k2 = mu*k1/((1+mu)^2);
c2 = 2*m2*sqrt(k2/m2)*sqrt(3*mu/(8*(1+mu))),
Ai = [m1*s^2+(c1+c2)*s+(k1+k2), -(c2*s+k2);
    -(c2*s+k2), m2*s^2+c2*s+k2];
Bi = [1;0];
H2 = inv(Ai)*Bi;     H2 = minreal(H2);
H21 = H2(1);         H22 = H2(2);
% ------- Amplitudengänge
fmin = f01/10;              fmax = f01*10;
e1 = round(log10(fmin));    e2 = round(log10(fmax));
f = logspace(e1,e2,5000);   w = 2*pi*f;
figure(1),   clf;   %##############
subplot(211),
[M11, phi11, w11] = bode(H11,w);
[M21, phi21, w21] = bode(H21,w); % Normaler Tilger
semilogx(w11/(2*pi), 20*log10(squeeze(M11)));
hold on;
semilogx(w11/(2*pi), 20*log10(squeeze(M21)),'r');
hold off;
title(['Frequenzgang Y1/Fex mit R =', num2str(R),...
    ' Ohm und Y1/Fex mit normaler Tilgung']);
grid on;   xlabel('Hz');
legend('Y1/Fex mit R', 'Y1/Fex');
subplot(212),
[M12, phi12, w12] = bode(H12,w);
[M22, phi22, w22] = bode(H22,w);
semilogx(w12/(2*pi), 20*log10(squeeze(M12)));
hold on;
semilogx(w12/(2*pi), 20*log10(squeeze(M22)),'r');
hold off
title(['Frequenzgang Y2/Fex mit R =', num2str(R),...
    ' Ohm und Y2/Fex mit normaler Tilgung']);
grid on;   xlabel('Hz');
legend('Y2/Fex mit R', 'Y2/Fex');

% -------- Schätzung der elektrischen Leistung
% wenn die Anregung mit einer Frequenz w01 stattfindet
% Frequenzgänge für w = j*w01
H11_w01 = polyval(b11, j*w01)/polyval(a11, j*w01);
H12_w01 = polyval(b12, j*w01)/polyval(a12, j*w01);
H13_w01 = polyval(b13, j*w01)/polyval(a13, j*w01);

Fex_ampl = 1,
em_ampl = abs(kv*H13_w01)*Fex_ampl,
i_ampl = em_ampl/R,
P = em_ampl*i_ampl/2,

% -------- Amplituden der Bewegung
y1_ampl = abs(H11_w01)*Fex_ampl,
y2_ampl = abs(H12_w01)*Fex_ampl,

% -------- Amplitude der Fem
Femf_ampl = kf*abs(H13_w01)*Fex_ampl,

% ----Zustandsmodell für C = 0 und i als Zustandsvariable
Az = [0 1 0 0 0;
    -(k1+k2)/m1,-c1/m1,k2/m1,0,kf/m1;
      0 0 0 1 0;
      k2/m2,0,-k2/m2,0,-kf/m2;
      0,-kv/L,0,kv/L,-R/L];
Bz = [0 1/m1 0 0 0]';
Cz = [1 0 0 0 0];             Dz = 0;
eigenwerte = eig(Az),



