% Skript rad_aufhaeng_3.m, in dem die Radaufhängung
% eines Fahrzeugs simuliert wird
% Es wird auch die partikuläre Lösung ermittelt
% Arbeitet mit Simulink-Modell rad_aufhaeng3.mdl

clear;
s = tf('s');
% ------- Parameter des Systems
ms = 300;   c = 12000;     k = 50000;
ma = 30;    ct = 10000;    kt = 100000;
f10 = sqrt(k/ms)/(2*pi),   % Eigenfrequenzen ohne Dämpfung
f20 = sqrt(kt/ma)/(2*pi),
% ------- Matrizen des Systems
M = [ms 0;0 ma];   C = [c -c;-c ct+c];   K = [k -k; -k kt+k];
% ------- Poly-Eigenwerte
[X,e] = polyeig(K, C, M);
[e, k] = sort(e),
X = X(:,k);
f1 = abs(imag(e(1)))/(2*pi),    f2 = abs(imag(e(3)))/(2*pi),

% ------- Homogene Lösung
dt = 0.01;       Tfinal = 1;
t = 0:dt:Tfinal;
xv0 = [0.2; 0; 0; 0]; % Anfangsbedingungen (Lagen und Geschwindigkeiten)
p = inv([X(:,1),X(:,2),X(:,3),X(:,4);...
    X(:,1)*e(1),X(:,2)*e(2),X(:,3)*e(3),X(:,4)*e(4)])*xv0,

xh_temp = p(1)*X(:,1)*exp(e(1)*t) + p(2)*X(:,2)*exp(e(2)*t)+...
    p(3)*X(:,3)*exp(e(3)*t) + p(4)*X(:,4)*exp(e(4)*t);
xh = real(xh_temp);     % Die eventuellen kleinen Imaginärteile
                        % wegen numerischen Fehlern entfernen
figure(1);    clf;
plot(t, xh);
title('Homogene Lösung für die Lagen der zwei Massen');
grid on;  xlabel('Zeit in s');     legend('x1', 'x2');

% ------- Partikuläre Lösung über Übertragungsfunktion
D = [0;ct];     B = [0;kt];
H = inv(s^2*M + s*C + K)*[s*D + B],
H1 = H(1),      H2 = H(2), % Transfer-Function Daten

[b, a] = tfdata(H);
b1 = b{1};    a1 = a{1};   % Koeffizienten der Übertragungs-
                           % funktion von y zu x1 
b2 = b{2};    a2 = a{2};   % Koeffizienten der Übertragungs-
                           % funktion von y zu x2
% -------- Aufruf der Simulation                           
dt = 0.01;    Tfinal = 1; 
my_options = simset('MaxStep',dt);
sim('rad_aufhaeng3',[0,Tfinal], my_options);                           

t = x.time;
x1 = x.signals.values(:,1);
x2 = x.signals.values(:,2);

figure(2);    clf;
plot(t, [x1, x2]);
title('Homogene plus partikuläre Lösung');
xlabel('Zeit in s');     grid on;
legend('x1(t)','x2(t)');

% -------- Frequenzgänge
fx = logspace(-1, 2); 
H1s = freqs(b1, a1, 2*pi*fx);
H2s = freqs(b2, a2, 2*pi*fx);
                           
figure(3);    clf;
semilogx(fx, 20*log10(abs(H1s)));
hold on;
semilogx(fx, 20*log10(abs(H2s)),'r');
title('Amplitudengänge')
hold off;    xlabel('Hz');   grid on; 
legend('y  zu  x1','y  zu  x2'); ylabel('dB')

                           
                           
                           
                           