% Skript uebertragung_mechanik_1.m, in dem die Übertragungsfunktionen eines
% mechanischen Systems mit einer Pulsanregung identifiziert wird

clear;
s = tf('s');
% -------- Parameter des Feder-Masse-System
m2 = 5;      k2 = 10;           c2 = 0.3;
% Optimale Tilgerparameter nach Hartog
m1 = 2;      mu = m1/m2;        k1 = mu*k2/((1+mu)^2);
w1 = sqrt(k1/m1);               c1 = 2*w1*m1*sqrt(3*mu/(8*(1+mu)^3));
% Ohne Tilger
%m1 = 2;      k1 = 0;            c1 = 0; 
% -------- Zustandsmodell
A = [-c1/m1, -k1/m1, c1/m1, k1/m1;
    1 0 0 0;
    c1/m2, k1/m2, -(c1+c2)/m2 -(k1+k2)/m2;
    0 0 1 0];
B = [0, 0, 1/m2, 0]'; 
C = zeros(4,4);    C(2,2) = 1;     C(4,4) = 1;  % Ausgangsvariablen y1, y2
D = zeros(4,1);
[eig_vekt, lambda] = eig(A);    % Eigenvektoren und Eigenwerte
lambda = diag(lambda);
f1 = imag(lambda(1))/(2*pi),
f2 = imag(lambda(3))/(2*pi),
fmax = max([f1, f2]);    Tmin = 1/fmax;

% Definition eines Zustandssystem (ss-System)
my_system = ss(A, B, C, D);

% Übertragungsfunktionen für die Ausgänge
H = C*inv(s*eye(4,4)-A)*B + D,
% -------- Komplexer Frequenzgang
f = logspace(-1, 1, 500);    w = 2*pi*f;
[mag, phi] = bode(H,w);

betrag1 = squeeze(mag(2,1,:));        % Amplitudengang für y1 
phase1 = squeeze(phi(2,1,:));         % Phasengang

betrag2 = squeeze(mag(4,1,:));        % Amplitudengang für y2
phase2  = squeeze(phi(4,1,:));         % Phasengang

figure(1);    clf;
subplot(411), semilogx(f, 20*log10(betrag1));
title('Amplitudengang von u(t) bis y1(t)');
xlabel('Hz');      ylabel('dB');     grid on,
subplot(412), semilogx(f, phase1);
title('Phasengang von u(t) bis y1(t)');
xlabel('Hz');      ylabel('Grad');     grid on,

subplot(413), semilogx(f, 20*log10(betrag2));
title('Amplitudengang von u(t) bis y2(t)');
xlabel('Hz');      ylabel('dB');     grid on,
subplot(414), semilogx(f, phase2);
title('Phasengang von u(t) bis y2(t)');
xlabel('Hz');      ylabel('Grad');     grid on,




