% Skript unwucht_rueckw_1.m, in dem die Schwingungen durch Unwucht
% eines Feder-Masse-System mit Motorantrieb untersucht werden
% Arbeitet mit Modell unwucht_rueckw1.slx
clear;
s = tf('s');
% ------- Parameter des Systems
m = 100;     m0 = 1;
kf = 2000;   c = 4;    
J = 10;      cj = 0.1;
m0e = 1;
U = 50;      km = 1;
g = 9.89;
L = 1;       R = 0.1;
% ------- Parameter des Systems
m = 100;     m0 = 1;
kf = 2000;   c = 4;    
J = 10;      cj = 0.1;
m0e = 1;
U = 50;      km = 10;     % Schwacher Motor
g = 9.89;
L = 1;       R = 0.1;
% ------- Eigene Frequenz ohne Daempfung
omega_0 = sqrt(kf/(m+m0)),
f_0 = omega_0/(2*pi),

% ------- Aufruf der Simulation
dt = 0.002;     Tfinal = 200;
t = 0:dt:Tfinal;
my_options = simset('MaxStep', dt);
sim('unwucht_rueckw1', [t], my_options);

%return
t = xy.time;
nt = length(t),
x = xy.signals.values(:,1);
omega = xy.signals.values(:,2); 
i = xy.signals.values(:,3);
theta = xy.signals.values(:,4);


figure(1);    clf;
subplot(411), plot(t, x);
title(['Bewegung der Hauptmasse']);
xlabel('Zeit in s');    grid on;    

subplot(412), plot(t, omega);
title('Kreisgeschwindigkeit omega');
xlabel('Zeit in s');    grid on;

subplot(413), plot(t, i);
title('Strom des Gleichstrommotors');
xlabel('Zeit in s');    grid on;

subplot(414), plot(t, theta);
title('Winkel der Drehachse');
xlabel('Zeit in s');    grid on;

% ------- Spektrum der Bewegung der Hauptmasse
h = spectrum.welch;          % Create a Welch spectral estimator. 
h.SegmentLength = 512*128;
h.SegmentLength = 512*64;
Fs = 1/dt;

nsta = fix(nt/2);            % Ohne Einschwingen
x = x(nsta:end);
omega = omega(nsta:end);
i = i(nsta:end);


Hpsd = psd(h,x,'Fs',Fs);             % Calculate the PSD 
Hpsd1 = psd(h,omega,'Fs',Fs);        % Calculate the PSD 
Hpsd2 = psd(h,i,'Fs',Fs);            % Calculate the PSD 

figure(2);     clf;
nd = 1000;
subplot(211), plot(Hpsd.Frequencies(1:nd), 10*log10(Hpsd.Data(1:nd)));                           
title('Spektrale Leistungsdichte der Bewegung der Hauptmasse x');
xlabel('Hz');    grid on;    ylabel('dB Watt/Hz');  

subplot(212), plot(Hpsd.Frequencies(1:nd), 10*log10(Hpsd1.Data(1:nd)));                           
title('Spektrale Leistungsdichte der Kreisgeschwindigkeit omega');
xlabel('Hz');    grid on;    ylabel('dB Watt/Hz');    

fdreh = mean(omega)/(2*pi),








