% Skript unwucht_synchro1.m,in dem das Modell
% unwucht_synchro_1.mdl initialisiert und aufgerufen wird
clear;
% ------- Parameter des Systems
M = 5;
M = 1000;
b = 1.25;          c = 20*pi*pi;
f0 = sqrt(c/M)/(2*pi),
% ------- Motore und Unwuchten
J1 = 1e-5;           
J2 = J1;
m1 = 0.1;            m2 = 0.1;
epsilon1 = 0.01;     epsilon2 = epsilon1;
k1 = 0.5e-4;         L10 = 25e-4;         
k2 = k1;             L20 = 24.5e-4;
%k2 = 0;              L20 = 0;  % Motor 2 ausgeschaltet
% ------- Aufruf der Simulation
Tfinal = 100;
dt = 0.001
t = 0:dt:Tfinal;
phi10 = 0;      dphi10 = 20;    % Anfangsbedibgungen
phi20 = 0;      dphi20 = 0;
%sim('unwucht_synchro_1', t);
sim('unwucht_synchro_2', t);

omega1 = y.signals.values(:,1);
omega2 = y.signals.values(:,2);
x = y.signals.values(:,3);
phi1 = phi.signals.values(:,1);
phi2 = phi.signals.values(:,2);
t = y.time;

figure(1),    clf;
subplot(211), plot(t, [omega1, omega2]);
title('Kreisgeschwindigkeiten');
xlabel('Zeit in s');   grid on;

subplot(212), plot(t, x);
title('Koordinate x des Schwingers');
xlabel('Zeit in s');   grid on;

figure(2),    clf;
plot(t, mod(phi1 - phi2,2*pi));
title('Pasendifferenz Modulo 2*pi');
xlabel('Zeit in s');    grid on;


% ------- Spektrale Analyse
nx = length(x);
nf = fix(nx/2);    % Einschwingen entfernen

O1 = fft(detrend(omega1(nf:end)));
O2 = fft(detrend(omega2(nf:end)));
X = fft(detrend(x(nf:end)));
nfft = length(X);

figure(3),    clf;
nd = fix(nfft/20);
subplot(311), plot((0:nd-1)/(dt*nfft), 20*log10(abs(O1(1:nd))/nfft));
title('Amplitudenspektrum der Schwankungen der Kreisgeschwindigkeit 1');
xlabel('Hz');    ylabel('dB');     grid on;
La = axis;   axis([La(1:2), -150, La(4)]);
subplot(312), plot((0:nd-1)/(dt*nfft), 20*log10(abs(O2(1:nd))/nfft));
title('Amplitudenspektrum der Schwankungen der Kreisgeschwindigkeit 2');
xlabel('Hz');    ylabel('dB');     grid on;
La = axis;   axis([La(1:2), -150, La(4)]);
subplot(313), plot((0:nd-1)/(dt*nfft), 20*log10(abs(X(1:nd))/nfft));
title('Amplitudenspektrum der Schwankungen der Koordinate x der Hauptmasse');
xlabel('Hz');    ylabel('dB');     grid on;
La = axis;   axis([La(1:2), -200, La(4)]);


