close all
clear all
%% Parameterdefinition
g=9.81;                     %[m/s^2]
m1=1;                       %[kg]
m2=5;                       %[kg]
m3=0.5;                       %[kg]
k=300;                     %[N/m]
k_phi=100;                 %[Nm]
c=0.4;                      %[kg/s]
c_phi=0.01;                  %[kgm^2/s]
l1=0.2;                     %[m]
M01=0.1;
F1=1;

%% Bestimmung der Eigenfrequenzen

%Steigigkeits- u. Massenmatrix

K=[k 0;0 k_phi+g*l1*(0.5*m1+m3)];
C=[c 0;0 c_phi];
M=[m1+m2+m3 0;0 l1^2*(1/3*m1+m3)];

% Eigenwerte (EW) und Eigenvektoren (EV) berechnen

[EV,EW]=eig(K,M);

% Eigenkreisfrequenzen

om0_1= sqrt(EW(1,1));
om0_phi=sqrt(EW(2,2));

% Eigenfrequenzen

f0_1 = om0_1 / (2*pi);
f0_phi = om0_phi / (2*pi);

% Zustandsraum

A=[-c/(m1+m2+m3) 0 -k/(m1+m2+m3) 0
    0 -c_phi/(l1^2*(1/3*m1+m3)) 0 -(k_phi+g*l1*(0.5*m1+m3))/(l1^2*(1/3*m1+m3))
     1 0 0 0
     0 1 0 0];
B=[1/(m1+m2+m3) 0;0 1/(1/3*m1+m3);0 0;0 0];
C=[0 0 1 0;0 0 0 1];
D=[0 0;0 0];

sys1=ss(A,B,C,D);

[M1,Phase,om]=bode(sys1,logspace(-1,2,2000));

bode(sys1)
hold on
f0=om(:)/(2*pi);

figure(1);
plot(f0(:),M1(:),'k--');
xlabel('Hz'); ylabel('H(\omega)')
title('frequency response function')
grid on
xlim([0 10])
ylim([0 0.3])