%   Script NULLST.M
%
%   Demonstration der Wirkung von Nullstellen
%   am Beispiel einer Anordnung mit drei Federn und
%   zwei Massen
%
% J. Lunze
% 23.2.1996
% für 2. Auflage: 14.12.1998
% für 3. Auflage: 3.4.2001
% für 4. Auflabe: 15.9.2003
% für 5. Auflage: 18.4.2005
% für 7. Auflage: 28.4.2008
% für 8. Auflage: 24.1.2010
% für 9. Auflage: 13.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 23.03.2020
%
%   
clear
close all
BspSchwinger;     %  Modell des Feder-Masse-Dämpfer-Systems
                  %  Schwinger=ss(A, b, c, d)
%
% Nullstellen des Systems:
tzero(Schwinger)
% Eigenwerte des Systems;
eig(A)
%
Tend = 70;
T = [0:0.1:Tend];   % Zeitskala mit 701 Punkten


%%                  
%
%   Feder-Masse-Schwinger bei u(t) = exp(0,0344t) sin(1,206t)
%   d.h. Erregung "in der Nullstelle"
%
figure(1);
subplot(6,1,[1,2]);
for k=0:700,
t = k/10;
U1(k+1,1)=exp(0.0344*t) * sin(1.2056*t);
   end;
plot(T, U1, 'b', [-5 70], [0 0], 'k:');
hold on;
plot([-5 0 0],[0 0 U1(1,1)], 'b');
xleer;
axis([-5, 70, -10, 10]);
yticks([-10 0 10]);
%latextitle('Eingangsgr\"oße $$u(t) = exp(0,0344 t) sin(1,2056 t)$$');
latexylabel('$$u$$ in N');
rotateY;
%
subplot(6,1,[3,4]);
[Y1,T,X1]=lsim(Schwinger, U1, T);
plot(T, Y1, 'b');
%latextitle('Systemverhalten bei Erregung mit $$u(t) = exp(0,0344 t) sin(1,2056 t)$$')
%xlabel('Zeit');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 Y1(1,1)], 'b');
latexylabel('$$y$$ in cm');
rotateY;
axis([-5 70 -2 2]);
yticks([-2 0 2]);
xleer;
hold off;
%
subplot(6,1,[5,6]);
plot(T, X1(:,2), 'b');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 X1(1,2)], 'b');
latexxlabel('$$t$$ in s');
latexylabel('$$x_2$$ in cm');
rotateY;
axis([-5 70 -20 20]);
xticks([0 20 40 60 80]);
yticks([-20 0 20]);
hold off;


%% Vergleich mit rein sinusförmiger Erregung
%
figure(2);
for k=0:700
t = k/10;
U2(k+1,1)=sin(1.2056*t);
end
subplot(6,1,[1,2]);
plot(T, U2, 'b', [-5 70], [0 0], 'k:');
hold on
plot([-5 0 0],[0 0 U2(1,1)], 'b')
axis([-5, 70, -1.5, 1.5]);
yticks([-1 0 1]);
%latextitle('Eingangsgr\"oße $$u(t) = sin(1,2056 t)$$');
xleer;
latexylabel('$$u$$ in N');
rotateY;

subplot(6,1,[3,4]);
Y2=lsim(Schwinger, U2, T);
plot(T, Y2, 'b');
%latextitle('Systemverhalten bei Erregung mit $$u(t) = sin(1,2056 t)$$')
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 Y2(1,1)], 'b' );
latexxlabel('$$t$$ in s');
latexylabel('$$y$$ in cm');
rotateY;
axis([-5 70 -1.5 1.5]);
xticks([0 20 40 60 80]);
yticks([-1 0 1]);
hold off;
 

%%   Anfangszustand, so dass die Eigenbewegung das Übergangsverhalten
%    kompensiert
%    
figure(3)
subplot(3,1,1);
s01=0.0344+1.2056j;
s02=0.0344-1.2056j;
u1=1/2/j;
u2=-u1;
%    Anfangszustand, für den die Eigenbewegung das Übergangsverhalten
%    kompensiert:
x0=inv(s01*eye(5,5)-A)*b*u1+inv(s02*eye(5,5)-A)*b*u2;
plot(T, U1, 'b', [-5 70], [0 0], 'k:');
hold on;
plot([-5 0 0],[0 0 U1(1,1)], 'b');
xleer;
axis([-5, 70, -10, 10]);
yticks([-10 0 10]);
latextitle('Systemverhalten bei Nulldynamik')
latexylabel('$$u$$ in N');
rotateY;

subplot(3,1,2);
[Y10,T,X10]=lsim(Schwinger, U1, T, x0);
plot(T, Y10, 'b');
%latextitle('Systemverhalten bei Erregung mit $$u(t) = exp(0,0344 t) sin(1,2056 t)$$')
%xlabel('Zeit');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 Y10(1,1)], 'b');
latexylabel('$$y$$ in cm');
rotateY;
axis([-5 70 -2 2]);
yticks([-2 0 2]);
xleer;
hold off;

subplot(3,1,3);
plot(T, X10(:,2), 'b');
%latextitle('Systemverhalten bei Erregung mit $$u(t) = exp(0,0344 t) sin(1,2056 t)$$')
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 X10(1,2)], 'b');
latexxlabel('$$t$$ in s');
latexylabel('$$x_2$$ in cm');
rotateY;
axis([-5 70 -20 20]);
xticks([0 20 40 60 80]);
yticks([-20 0 20]);
hold off;



%%   y(t)=0 bei beliebigem Anfangszustand xN0
%
figure(4)
%  Transformation des Systems in E/A-Normalform
SS=ctrb(A, b);
SSinv=inv(SS);
sRT=SSinv(n,:);
TEinv=[c; c*A; sRT; sRT*A; sRT*A*A];
TE=inv(TEinv);
AE=TEinv*A*TE;
bE=TEinv*b;
cE=c*TE;
EANormalform=ss(AE, bE, cE, d);
T=0:0.1:70;
ANN=AE(3:5, 3:5);
bN=bE(3:5);
bq=c*A*b;
cN=1/bq*AE(2, 3:5);  % Ausgabe von yN=cN xN
dN=0;
uGenerator=ss(ANN, bN, cN, dN);
% Transformation in kanonische Normalform, um Anfangsbedingungen
% für spezielle Signale auswählen zu können
%  TMatrix = Transformationsmatrix, so dass xtilde=TMatrix*xN
[uGeneratorKanon,TMatrix]=canon(uGenerator, 'modal');
%  2. und 3. Zustandsvariable bestimmen die aufklingende Sinusfunktion
%  xN0tilde=-4*[0; 1; 1];  liefert ungefähr die Eingangsgröße U1
xN0tilde=250000000*[1;0;0];
xN0=inv(TMatrix)*xN0tilde;
YN=initial(uGeneratorKanon, xN0tilde, T);
UE=-YN;
subplot(3,1,1)
plot(T, UE, 'b',  [-5 70], [0 0], 'k:');
hold on;
plot([-5 0 0],[0 0 UE(1,1)]);
xleer;
axis([-5, 70, -1, 1]);
yticks([-1 1]);
%latextitle('Eingangsgr\"oße $$u(t) = exp(0,0344 t) sin(1,2056 t)$$');
latexylabel('$$u$$ in N');
rotateY;
latextitle('Systemverhalten bei Nulldynamik')
hold off;

subplot(3,1,2);
xE0=[0;0;xN0];
[YE, T, XE]=lsim(EANormalform, UE, T, xE0);
plot(T, YE, 'b');
%latextitle('Systemverhalten bei Erregung mit $$u(t) = exp(0,0344 t) sin(1,2056 t)$$')
%xlabel('Zeit');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 YE(1,1)], 'b');
latexylabel('$$y$$ in cm');
rotateY;
axis([-5 70 -1 1]);
yticks([-1 1]);
xleer;
hold off;

subplot(3,1,3);
cE2=[0 1 0 0 0]*TE;    % Ausgabe von x_2 aus Originalmodell
X2=XE*cE2';
plot(T, X2, 'b');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 X2(1,1)], 'b');
latexxlabel('$$t$$ in s');
latexylabel('$$x_2$$ in cm');
rotateY;
axis([-5 70 -1 1]);
xticks([0 20 40 60 80]);
yticks([-1 1]);
hold off;
 

%%    Dasselbe für anderen Anfangszustand xN0
% 
figure(5)

%  2. und 3. Zustandsvariable bestimmen die aufklingende Sinusfunktion
%  xN0tilde=-4*[0; 1; 1];  liefert ungefähr die Eingangsgröße U1
%xN0tilde=[2500; -4; 4];
xN0tilde=100000*[2500; -4; 4];
xN0=inv(TMatrix)*xN0tilde;
YN=initial(uGeneratorKanon, xN0tilde, T);
UE=-YN;
%
subplot(6,1,[1,2])
plot(T, UE, 'b',  [-5 70], [0 0], 'k:');
hold on;
plot([-5 0 0],[0 0 UE(1,1)], 'b');
xleer;
axis([-5, 70, -1, 1]);
yticks([-1 1]);
latexylabel('$$u$$ in N');
rotateY;
latextitle('Systemverhalten bei Nulldynamik')
hold off;
%
subplot(6,1,[3,4]);
xE0=[0;0;xN0];
[YE, T, XE]=lsim(EANormalform, UE, T, xE0);
plot(T, YE, 'b');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 YE(1,1)], 'b');
latexylabel('$$y$$ in cm');
rotateY;
axis([-5 70 -1 1]);
yticks([-1 1]);
xleer;
hold off;
%
subplot(6,1,[5,6]);
X2=XE*cE2';
plot(T, X2, 'b');
hold on
plot([-5 70], [0 0], 'k:', [-5 0 0],[0 0 X2(1,1)], 'b');
latexxlabel('$$t$$ in s');
latexylabel('$$x_2$$ in cm');
rotateY;
axis([-5 70 -1 1]);
xticks([0 20 40 60 80]);
yticks([-1 1]);
hold off;


%% Herstellung der Bilddateien *.eps
%
Epsbildc15('nullst');
