close all
clear variables

% symbolische Auswertung
%syms Ja Jp Js Jl;                                          % Massentraegheitsmomente
%syms d1 d2 da dt dl;                                       % Daempfungskonstanten
%syms c1 c2 ct cl;                                          % Federkonstanten
%syms Mp Ms;                                                % Schnittmomente an Momentaufnehmern

% Massenträgheiten in kg m²
Jp = 0.3;
Js = 0.2;
Jl = 0.9;

% Daempfungskonstanten in kg m² / s
da = 0.07;
dt = 0.02;
dw = 0;
dl = 0.07;

% Federkonstanten in N m / rad
cw = 10000;
ct = 250;

% Definitionsbereich der Plots
w = {1, 1000};

% Abbildung der Uebertragungsfunktion in Abhaengigkeit der
% Torsionssteifigkeit der Welle
while cw <= 500
    
% Matrizen des dynamischen Systems
v = [Jl Js Jp];                                             % Diagonale der Massenmatrix
M = diag(v);                                                % Massenmatrix
D = [dl+dw -dw 0;                                           % Daempfungsmatrix
    -dw dw+dt -dt;
    0 -dt dt+da];         
K = [cw -cw 0;                                              % Steifigkeitsmatrix
    -cw cw+ct -ct;
    0 -ct ct];            
F = [1 0 0;                                                 % Eingangsmatrix
    0 -1 0;
    0 1 -1];                 

% Zustandsraum
A = [zeros(3,3) eye(3,3);                                   % Systemmatrix
    -M\K -M\D];             
B = [zeros(3,3);                                            % Eingangsmatrix
    M\F];
C = [ct -ct 0 dt -dt 0;                                     % Ausgangsmatrix
    0 0 0 1 0 0;
    0 cw -cw 0 dw -dw;
    0 0 0 0 0 1];
d = [0 -1 0;                                                % Durchgangsmatrix
    0 0 0;
    0 0 0;
    0 0 0];                                             

sys = ss(A,B,C,d);                                          % Zustandsraumdarstellung
G = tf(sys);                                                % Uebertragungsfunktion

    bode(G(3,3),w);
    cw = cw+500;
    hold on
end
hold off

size(ss(A,B,C,d));                                          % Kontrolle