%   Script Syst2ord.m
%
%   Verhalten von Systemen zweiter Ordnung
%
% J. Lunze
% 9.10.95
% für 3. Auflage: 10.2.2001
% für 4. Auflage: 30.9.2003
% für 5. Auflage: 19.4.2005
% für 6. Auflage: 14.12.2006
% für 7. Auflage: 19.6.2008
% für 8. Auflage: 5.1.2010
% für 9. Auflage: 24.5.2011
% für 9. Auflage: 16.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 23.03.2020
%
%   
clear
close all
A = [-0.8 2;
     -2 -0.8];
b = [1;
     1];
cT = [1 1];
d = 0.3;
System1=ss(A, b, cT, d);
x0=[1, 1];

%%  Übergangsfunktion  
%   sprungfähigen Systems zweiter Ordnung
%   
figure(1);
subplot(2,1,1);
[Y, T, X]=step(System1);
plot(T, Y, 'b', [-1 0 0],[0 0 d],'b',[0 0],[-0.5 1.5],'k:', [-1 8],[0 0],'k:');
xleer;
yleer;
latexxlabel('$$t$$');
latexylabel('$$h$$');
rotateY;
latextitle('\"Ubergangsfunktion eines Systems zweiter Ordnung');
latextext(T(1,1)-0.5, Y(1,1), '$$d$$');
latextext(T(1,1),  Y(1,1),'-');
xticks([0]);
yticks([0]);
[m, n]=size(T);
Tmax = max(m, n);
axis([-1 T(Tmax) -0.2 1.3]);
latextext(T(Tmax)+0.2, Y(Tmax), '$$k_{\rm s}$$')
 

%%   Gewichtsfunktion
%   (Sprungfähigkeit in Gewichtsfunktion nicht sichtbar)
%
figure(2);
subplot(2,1,1);
[Y, T, X]=impulse(System1);
plot(T, Y, 'b', [-1 0 0],[0 0 cT*b], 'b', [0 0],[-1 2.5],'k:',[-1 7],[0 0],'k:');
xleer;
yleer;
%latexxlabel('$$t$$');
latexylabel('$$g$$');
rotateY;
latextitle('Gewichtsfunktion eines Systems zweiter Ordnung');
xticks([0]);
yticks([0]);
axis([-0.7 7 -0.9 2.3]);
latextext(5.6, 1.8, '$$r=1$$');
%   System mit g(0)=0
PT2Glied=ss([-0.1 0; 1 -0.3], [1; 0], [0 2*0.3*0.1], 0);
subplot(2,1,2);
[Y2, T2, X2]=impulse(PT2Glied);
plot(T2, Y2, 'b', [-7 0],[0 0], [0 0],[-0.02 0.2],'k:',[-7 70],[0 0],'k:');
xleer;
yleer;
latexxlabel('$$t$$');
latexylabel('$$g$$');
rotateY;
xticks([0]);
yticks([0]);
axis([-6.5 65 -0.01 0.14]);
latextext(52, 0.11, '$$r=2$$');
 

%%  Übergangsfunktionen von P-Gliedern
%   
figure(3);
subplot(2,1,1)
T3=0:0.05:50;
PGlied=ss([0], [0], [0], 2);
PT1Glied=ss([-0.3], [1], [2*0.3], 0);
Y3=step(PT2Glied, T3);
Y1=step(PGlied, T3);
Y2=step(PT1Glied, T3);
plot(T3, Y1, 'b', T3, Y2, 'b', T3, Y3, 'b', [-10 50],[0 0],'k:',...
          [-10 50],[2 2],'k:', [-10 0 0],[0 0 2], 'b',[0 0],[-0.1 2.3],'k:');
xleer;
yleer;
latexxlabel('$$t$$');
latexylabel('$$h$$');
rotateY;
latextitle('\"Ubergangsfunktionen');
%latextext(1, 2.15, 'P-Glied');
%latextext(8, 1.7, 'PT_1-Glied');
%latextext(20, 1.2, 'PT_2-Glied ($$d>1$$)')
xticks([0]);
yticks([0]);
latextext(-15, 2, '$$k_{\rm s}$$');
axis([-10 50 -0.1 2.3]);


%%   Stochastisch gestörte Übergangsfunktion eines PT1-Gliedes
%   
figure(4);
dt=0.2;
Te=20;
T=[0:dt:Te];
System=ss([-0.3], [1], [4.3*0.3], 0);
[Y2, X]=step(System, T);
F=rand(size(Y2));
Y2 = Y2 + F;
Tt=5;
kneu=Tt/dt;
Y2t=zeros((Te+Tt)/dt,1);
for k=1:length(T)
    Y2t(k+kneu,1)=Y2(k);
    Tneu(k+kneu)=T(k);
end;
Tneu(1:kneu)=[-Tt:dt:-dt];
Y2t(1:kneu,1)=rand(kneu,1);
plot(Tneu, Y2t, 'b');
latexxlabel('$$t$$ in min');
latexylabel('$$h$$ in m');
rotateY;
latextitle('Experimentell ermittelte \"Ubergangsfunktion')
axis([-5 20, 0 5.5]);
xticks([-5 0 5 10 15 20]);
yticks([0 1 2 3 4 5 6]);


%%    Ausgangsgröße eines IT1-Gliedes
%   
figure(5);
A=[-3 0;
     1 0];
b=[3;
   0];
c=[0 1];
d=0;
System1=ss(A, b, c, d);
System2=tf(29, [1 1.5 29]);
System3=ss(-3, 6, 1, 0);
T=0:0.02:10;
for k=1:49
   U(1,k)=0;
end;
for k=50:200
   U(1,k)=1;
end;
for k=201:300
   U(1,k)=0;
end;
for k=301:400
   U(1,k)=-1;
end;
for k=401:501
   U(1,k)=0;
end;
[Y, T, X]=lsim(System1, U', T');
[Y1, T, X]=lsim(System2, U', T');
[Y2, T, X]=lsim(System3, U', T');
subplot(8, 1, [1,2]);
plot(T, U, 'b', [0 10], [0 0],'k:');
axis([0 10 -1.5 1.5]);
latexylabel('$$u$$');
rotateY;
yticks([-1 1]);
xleer;
latextitle('Verhalten einfacher \"Ubertragungsglieder');
%
subplot(8, 1, [3,4]);
plot(T, Y2, 'b', [0 10],[0 0],'k:');
latexylabel('$$y$$');
rotateY;
axis([0 10 -2.2 2.2]);
yticks([-2 2]);
xleer;
%
subplot(8, 1, [5,6]);
plot(T, Y1, 'b', [0 10],[0 0],'k:');
latexylabel('$$y$$');
rotateY;
axis([0 10 -2 2]);
yticks([-2 2]);
xleer;
%
subplot(8, 1, [7,8]);
plot(T, Y, 'b', [0 10],[0 0],'k:');
latexylabel('$$y$$');
rotateY;
set(gca, 'YTick', [-1 3]);
rotateY;
latexxlabel('$$t$$');
axis([0 10 -0.5 4]);
yticks([0 3]);
xticks([0 2 4 6 8 10]);


%%   Beispiele für Systeme zweiter Ordnung
%   
figure(6)
A=[ 0  1  0   0 0   0  0  0  0  0;
  -18 -6  0   0 0   0  0  0  0  0;
    0  0  0   1 0   0  0  0  0  0;
    0  0 -34 -6 0   0  0  0  0  0;
    0  0  0   0 -2  0  0  0  0  0;
    0  0  0   0 1 -2.3 0  0  0  0;
    0  0  0   0 0   0 -5  0  0  0;
    0  0  0   0 0   0  1 -6  0  0;
    0  0  0   0 0   0  0  0  2  0;
    0  0  0   0 0   0  0  0  1  3];
B=[0 0  0  0  0;
   1 0  0  0  0;
   0 0  0  0  0;
   0 1  0  0  0;
   0 0 6.9 0  0;
   0 0  0  0  0;
   0 0  0 30  0;
   0 0  0  0  0;
   0 0  0  0 0.05;
   0 0  0  0  0 ];
C=[18 0 0  0 0 0 0 0 0 0;
    0 0 34 0 0 0 0 0 0 0;
    0 0 0  0 0 1/1.5 0 0 0 0;
    0 0 0  0 0 0 0 1 0 0;
    0 0 0  0 0 0 0 0 0 1];
D=zeros(5,5);
subplot(2,3,1);
Sys=ss(A, B, C, D);
[P, Z]=pzmap(Sys);
for i=1:length(P)
    plot(real(P(:)), imag(P(:)), 'bx');
    hold on
end
for i=1:length(Z)
    plot(real(Z(:)), imag(Z(:)), 'bo');
end
plot([-10 10],[0 0],'k:', [0 0 ],[-10 10],'k:');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('PN-Bild');
axis([-7 5 -6 6]);
xticks([-5 0 5]);
yticks([-5 0 5]);
axis('square');
T=0:0.02:2;
[Y, T, X]=step(Sys, T);
for i=1:5
   subplot(2, 3, i+1);
   plot(T, Y(:,i,i), 'b');
   axis([0 2 0 1.2]);
   xticks([0 1 2]);
   yticks([0 0.5 1 1.5]);
end;
%
subplot(2, 3, 2);
latextitle('$$h(t)$$');
subplot(2, 3, 4);
latexxlabel('$$t$$');
latexylabel('$$h$$');
rotateY;
subplot(2,3,5);
latexxlabel('$$t$$');
subplot(2,3,6);
latexxlabel('$$t$$');



%%   Übergangsfunktion eines Systems zweiter Ordnung
%   für d=0.5 und d=0.2
%   
figure(7);
A = [-0.8 2;
     -2 -0.8];
b = [1;
     1];
cT = [1 1];
d = 0.3;
x0=[1, 1];
T=0:0.1:40;
subplot(6,1,[1,2]);
System=tf([1], [1 0.5 1]);
[Y, X]=step(System, T);
plot(T, Y, 'b', [-6 0],[0 0], 'b',[0 0],[-0.5 2],'k:', [-2 30],[0 0],'k:');
xleer;
yleer;
latexylabel('$$h$$');
rotateY;
latextitle('\"Ubergangsfunktion eines Systems zweiter Ordnung');
xticks([0]);
yticks([0]);
latextext(15, 0.5, '$$d=0,\!5$$');
axis([-4 25 -0.3 1.8]);
%
subplot(6,1,[3,4]);
System2=tf([1], [1 0.2 1]);
[Y, X]=step(System2, T);
plot(T, Y, 'b', [-6 0],[0 0], 'b',[0 0],[-0.5 2],'k:', [-6 40],[0 0],'k:');
xleer;
yleer;
latexxlabel('$$t$$');
latexylabel('$$h$$');
rotateY;
yticks([0]);
xticks([0]);
latextext(15, 0.5, '$$d=0,\!2$$');
axis([-4 25 -0.3 1.8]);
 

%%   Übergangsfunktion von PT_2-Gliedern mit derselben Eigenfrequenz w_0=1
%   
figure(8)
T1=1;
T=0:0.1:30;
for d=[0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99]
    Nenner=[T1*T1, 2*d*T1, 1];
    System=tf(1, Nenner);
    Pole=roots(Nenner);
    subplot(2,3,3)
    plot(real(Pole(1,1)), imag(Pole(1,1)),'bx',  real(Pole(2,1)), imag(Pole(2,1)),'bx');
    hold on
    subplot(2,3,[1, 2]);
    Y=step(System, T);
    plot(T, Y, 'b');
    hold on
end
subplot(2,3,3)
k=0;
for phi=[90:5:270]
    k=k+1;
    x(k)=cos(phi/360*2*pi);
    y(k)=sin(phi/360*2*pi);
end
plot(x, y, 'k:', [-5 5],[0 0],'k:', [0 0],[-5 5],'k:');
axis('square');
latextitle(['Systeme mit der Zeitkonstanten $$T=$$',num2str(T1)]);
latexxlabel('Re');
rotateY;
axis([-1.1 1.1 -1.1 1.1]);
xticks([-1 0 1]);
yticks([-1 0 1]);
hold off
%
subplot(2,3,[1, 2]);
latexylabel('$$h$$');
rotateY;
latexxlabel('$$t$$');
latextext(0.1, 2.3, '$$d=$$');
latextext(0.1, 2, '0,01');
latextext(0.1, 1.7, '0,1');
latextext(0.1, 1.3, '0,3');
latextext(0.1, 1, '0,5');
axis([0 15 0 2.1]);
xticks([0 5 10 15]);
yticks([0 1 2]);
hold off
 

%% Übergangsfunktion von PT_2-Gliedern, deren Pole auf einer Parallelen
%  zur Imaginärachse liegen
%   
figure(9)
ddurchT=0.1;
subplot(2,3,3)
plot([-5 5],[0 0],'k:', [0 0],[-5 5],'k:', [-ddurchT -ddurchT],[10 -10],'k:');
hold on
T=0:0.1:30;
for d=[0.05,  0.2, 0.3, 0.4, 0.5]
    T1=d/ddurchT;
    Nenner=[T1*T1, 2*d*T1, 1];
    System=tf(1, Nenner);
    Pole=roots(Nenner);
    subplot(2,3,3)
    plot(real(Pole(1,1)), imag(Pole(1,1)),'bx',  real(Pole(2,1)), imag(Pole(2,1)),'bx');
    subplot(2,3,[1, 2]);
    Y=step(System, T);
    plot(T, Y, 'b');
    hold on
end
%   Einhüllende für die größte Verstärkung d
E1=1-exp(-ddurchT*T)/sqrt(1-d*d);
E2=1+exp(-ddurchT*T)/sqrt(1-d*d);
plot(T, E1,'b--', T, E2,'b--');
%  Beschriftung
subplot(2,3,3);
latextitle(['Systeme mit $$d/T= $$',num2str(ddurchT)]);
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-0.55 0.55 -0.55 0.55]);
xticks([-0.5 0 0.5]);
yticks([-0.5 0 0.5]);
axis('square');
hold off
%
subplot(2,3,[1, 2]);
latexylabel('$$h$$');
rotateY;
latexxlabel('$$t$$');
axis([0 30 -0.3 2.2]);
xticks([0 5 10 15 20 25 30]);
yticks([0 1 2]);
hold off
 

%%  Übergangsfunktion von PT_2-Gliedern mit unterschiedlichen Eigenfrequenz w_0=1
%   
figure(10)
d=0.4;
T=0:0.01:30;
for T1=[0.3, 0.5, 0.7, 0.9, 1]
    Nenner=[T1*T1, 2*d*T1, 1];
    System=tf(1, Nenner);
    Pole=roots(Nenner);
    subplot(2,3,3)
    plot(real(Pole(1,1)), imag(Pole(1,1)),'bx',  real(Pole(2,1)), imag(Pole(2,1)),'bx');
    hold on
    subplot(2,3,[1, 2]);
    Y=step(System, T);
    plot(T, Y, 'b');
    hold on
end
subplot(2,3,3)
plot([-5 5], [0 0],'k:', [0 0], [-5 5],'k:', [0 -5],[0 -90/40*5],'k:', [0 -5],[0 90/40*5],'k:');
axis('square');
latextitle(['Systeme mit der Zeitkonstanten $$T= $$',num2str(T1)]);
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-3.2 3.2 -3.2 3.2]);
xticks([-2 0 2]);
yticks([-2 0 2]);
hold off
%
subplot(2,3,[1, 2]);
latexylabel('$$h$$');
rotateY;
latexxlabel('$$t$$');
latextext(0.05, 1.1, '0,3');
latextext(1.7, 0.6, '$$T=1$$');
axis([0 6 0 1.3]);
xticks([0 2 4 6]);
yticks([0 1 2]);
hold off

%%  Ortskurven von PT_2-Gliedern mit derselben Eigenfrequenz w_0=1
%   
figure(11)
T1=1;
for d=[0.1, 0.3, 0.5, 0.7,  1]
    Nenner=[T1*T1, 2*d*T1, 1];
    System=tf(1, Nenner);
    [Real,Imag]=nyquist(System);
    plot(Real(:), Imag(:), 'b');
    hold on
end
plot([-3 4],[0 0],'k:', [0 0],[-6 6],'k:');
axis('square');
latextitle(['Systeme mit der Zeitkonstanten $$T= $$',num2str(T1)]);
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextext(0.2, -4.5, '$$d=0,\!1$$');
latextext(0.5, -2, '0,3');
latextext(0.2, -5.2, '$$\leftarrow$$');
latextext(0.5, -1.2, '0,5');
latextext(0.5, -0.2, '1');
plot(1, 0, 'bo');
latextext(0.9, 0.5, '$$k_{\rm s}=1$$');
axis([-3 3.5 -5.5 1]);
xticks([-2 0 2]);
yticks([-4 -2 0 2 4]);
hold off
 

%% Stationäres Verhalten eines PT2-Gliedes
%   
figure(12)
T1=1;
T2=0.5;
A=[-1/T1 0; 1/T2 -1/T2];
B=[1/T1; 0];
C=[0 1];
D=0;
PT2Glied=ss(A, B, C, D);
%   Eingangsgröße  u(t) = exp(mu*t)
mu=-1.5;
Tend=3.5;
T=0:0.01:Tend;
U=exp(mu*T);
x0=inv(mu*eye(2,2)-A)*B;
[Y, T]=lsim(PT2Glied, U, T, x0);
subplot(6,1,[1,2])
plot(T, U, 'b', [-0.3 0 0], [0 0 U(1,1)], 'b',[-0.3 Tend],[0 0],'k:');
latexylabel('$$u$$');
rotateY;
xleer;
axis([-0.3 Tend -0.1 1.1]);
yticks([0 1]);
latextitle('Station\"ares Verhalten eines PT2-Gliedes');
%
subplot(6,1,[3,4]);
plot(T, Y, 'b', [-0.3 0 0], [0 0 Y(1,1)], 'b',[-0.3 Tend],[0 0],'k:');
latexylabel('$$y_{\rm s}$$');
rotateY;
yticks([-8 -4 0]);
latexxlabel('$$t$$');
axis([-0.3 Tend -8.5 0.4]);
xticks([0 1 2 3 4]);
yticks([-8 -4 0]);
%  Berechnung der Amplitude des stationären Verhaltens
%  aus der Übertragungsfunktion
[z, n]=tfdata(PT2Glied,'v');
polyval(z,mu)/polyval(n,mu)


%%     Bilder auf .EPS- und .PS-Dateien
%
Epsbildc15('syst2or');
