%   Script NYQUKRIT.M
%
%   Beispiele zum Nyquistkriterium
%
% J. Lunze
% 20.1.1996
% für 3. Auflage: 15.5.2001
% für 4. Auflage: 3.10.2003
% für 5. Auflage: 13.6.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 19.6.2008
% für 9. Auflage: 16.2.2012
% für 10. Auflage: 27.1.2014
% für 11. Auflage: 18.4.2016
% für 12. Auflage: 24.07.2019
%
%   
clear
close all


%%
%  Beispiel 1:  Go = Kp(s+0.2)
%                    -------------------------------------
%                    (s^2 + 2s + 10)(s+4)(s^2 + 0.2s +0.1)
%   
figure(1);   % Bilder 1 und 10
z=[1 0.2];
n=[1 6.2 19.7 44.2 9.8 4];
System1=tf(z, n);
W=logspace(-2, 2, 1000);
%
subplot(4,4,[1 2 5 6]);
[Real, Imag]=nyquist(120*System1, W);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on;
plot([-1 -1],[-0.3 0.3],'k', [-1.5 -0.5],[0 0], 'k');
plot(Real(1), 0, 'bo');
axis([-3 17 -15 5]);
xticks([0 10]);
yticks([ -10  0 ]);
axis('square');
grid('on');
%latextitle('Offene Kette für $$k_P=120$$ und $$k_P=40$$');
latextext(2, 3, '$$\omega=0$$');
latextext(5.5, -12, '$$\leftarrow \omega$$');
latexxlabel('Re');
latexylabel('Im');
rotateY;
%
subplot(4,4,[3 4 7 8]);
[Real, Imag]=nyquist(40*System1);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on
plot([-1 -1],[-0.5 0.5],'k', [-1.5 -0.5],[0 0], 'k');
%plot(Real(1), 0, 'bo');
axis([-3 17 -15 5]);
xticks([0 10]);
yticks([-10 0 ]);
axis('square');
grid('on');
latextext(1, -6, '$$\leftarrow \omega$$');
latexxlabel('Re');
ylabel('');
hold off


%%  Beispiel 2:          Kp
%               Go = ----------    für  kp=4
%                    (s-1)(s+3)
%   
figure(2);
subplot(4,4,[1 2 5 6]);
z=[1];
n=[1 2 -3];
z=z*4;
System2=tf(z, n);
[Real, Imag]=nyquist(System2);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on;
plot([-1.05 -0.95], [0 0], 'k',[-1 -1], [0.05 -0.05],'k');
plot(Real(1), 0, 'bo');
grid('on');
latextitle('Ortskurve einer instabilen offenen Kette ');
latextext(-0.9, -0.5, '$$\omega \rightarrow$$');
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-1.4 0 -0.7 0.7]);
xticks([-1 0]);
yticks([-0.5 0.5]);
axis('square');

%%  Beispiel 2':          Kp(s+1)
%                  Go = ----------    für  kp=6
%                       (s-1)(s-3)
%   
subplot(4,4,[3 4 7 8]);
z=[1 1];
n=[1 -4 +3];
z=z*6;
System3=tf(z, n);
[Real, Imag]=nyquist(System3);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on
plot([-1.2 -0.8], [0 0], 'k', [-1 -1], [0.2 -0.2], 'k');
plot(Real(1), 0, 'bo');
grid('on');
latextitle('...mit zwei instabilen Polen ');
latextext(-0.5, 1.2, '$$\leftarrow \omega$$');
latexxlabel('Re');
ylabel('');
axis([-2.2 2.2 -2.2 2.2]);
xticks([-2 0 2]);
yticks([-2 0 2]);
axis('square');
hold off


%%  Beispiel 3:          Kp
%               Go = ----------    für  kp=4
%                    s(s+1)(s+3)  
figure(3);
z=[1];
n=[1 4 3 0];
z=z*4;
System4=tf(z, n);
subplot(1, 2, 1);
[Real, Imag]=nyquist(System4);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on;
axis([-2 4 -20 20]);
grid('on');
latextitle('Offene Kette mit I-Verhalten ');
%latextext(-1.02, 0, 'O');
axis([-2 1 -20 20]);
xticks([-2 -1 0 1 2]);
yticks([-20 -10 0 10 20]);
latextext(-0.5, 2.5, '$$\omega=-\infty$$');
latextext(-0.5, -2.5, '$$\omega=+\infty$$');
latextext(-1.6, 18, '$$\omega=-0$$');
latextext(-1.6, -18, '$$\omega=+0$$');
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold off
%
%   'Simulation' der Verbindungslinie für w->0
%   durch Verwendung eines Pols bei s0=-0.1 (anstelle bei 0)
%
z=[1];
s0 = 0.1;
n=[1 4+s0 3+4*s0 0+3*s0];
z=z*4;
System5=tf(z, n);
subplot(1, 2, 2);
[Real, Imag]=nyquist(z, n);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on
plot(Real(1), 0, 'bo');
axis([-3 14 -8.5 8.5]);
xticks([0 5 10]);
yticks([-5 0 5]);
axis('square');
grid('on');
latexxlabel('Re');
ylabel('');
latextext(4.5, -5, '$$\leftarrow \omega$$');



%%  Beispiel 4:          kp (1+0.4s)(1+0.2s)
%               Go = ---------------------------
%                    (s+0.1)(s+0.3)(s+0.5)(s+30)
%
%               für kp=5000  
figure(4);
z=5000*[0.08 0.6 1];
n=[1.0000   30.9000   27.2300    6.9150    0.4500];
System6=tf(z,n);
[Betrag, Phase, W]=bode(System6);
subplot(2,1,1)
semilogx(W, 20*log10(Betrag(:)), 'b');
latextitle('Offene Kette vierter Ordnung ');
xleer;
latexylabel('$$|G_0|$$ in dB');
rotateY;
grid('on');
axis([0.001 1000 -100 100]);
yticks([-100 -50 0 50 100]);
%
subplot(2,1,2);
semilogx(W, Phase(:), 'b');
latexxlabel('$$\omega$$ in rad/s');
latexylabel('$${\rm arg}(G_0)$$ in Grad');
rotateY;
axis([0.001 1000 -200 0]);
xticks([0.001 0.01 0.1 1 10 100 1000]);
yticks([-270 -180 -90 0]);
grid('on');


%%   Ortkurve desselben Beispiels
%
figure(5);
subplot(3,1,[1,2]);
[Real, Imag]=nyquist(System6);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
latextitle('Ortskurve eines Systems vierter Ordnung');
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-2000 12000 -10000 10000]);
xticks([0 5000 10000]);
yticks([-10000 0 10000]);
grid('on');


%%   Beispiel 5:          kp (s+0.48)(s+0.96)(s+1.44)(s+1.92)
%               Go = -----------------------------------------------
%                    (s+0.1)(s+0.2)(s+0.3)(s+0.4)(s+0.5)(s+20)(s+40)
%
%               für kp=1000  (d.h. ks= 8918)
%   
figure(6);
z=[7.0000   33.6000   56.4480   38.7072    8.9181];
n=[0.0010    0.0615    0.8908    1.2512    0.6935    0.1816    0.0220 0.0010];
System7=tf(z, n);
[Betrag, Phase, W]=bode(System7);
subplot(2,1,1);
semilogx(W, 20*log10(Betrag(:)), 'b');
hold on
semilogx([0.01 1000],[0 0],'k', [0.01 1000],[0 0],'b--');
xleer;
grid('on');
latextitle('Offene Kette siebenter Ordnung ');
latexylabel('$$|G_0|$$ in dB');
rotateY;
axis([0.01 1000 -100 100]);
yticks([-100 -50 0 50 100]);
hold off
%
subplot(2,1,2);
semilogx(W, Phase(:), 'b');
hold on
semilogx([0.01 1000],[-180 -180],'k', [0.01 1000],[-180 -180],'b--');
latexxlabel('$$\omega$$ in rad/s');
latexylabel('$$\phi_0$$ in Grad');
rotateY;
axis([0.01 1000 -290 0]);
xticks([0.01 1 100]);
yticks([-270 -180 -90 0]);
grid('on');
hold off


%%    Ortkurve desselben Beispiels
% 
figure(7);
subplot(4,4,[1 2 5 6]);
[Real, Imag]=nyquist(System7);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on
plot(Real(1), 0, 'bo');
plot([-3000 11000],[0 0],'k:', [0 0],[-7000 7000],'k:');
axis([-3000 11000 -7000 7000]);
xticks([0 5000 10000]);
yticks([-5000 0 5000]);
axis('square');
grid('on');
latexxlabel('Re');
latextitle('System siebenter Ordnung');
latexylabel('Im');
rotateY;
latextext(1500, -4000, '$$\leftarrow \omega$$');
%
subplot(4,4,[3 4 7 8]);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on
plot([-20 5],[0 0],'k:', [0 0],[-20 20],'k:');
plot([-1 -1],[-0.5 0.5],'k', [-1.5 -0.5],[0 0],'k');
grid('on');
axis([-15 1 -8 8]);
xticks([-15 -10 -5 0]);
yticks([-5 0 5]);
axis('square');
latextext(-9, -3,'$$\omega \rightarrow \infty$$');
latexxlabel('Re');
ylabel('');


%%   Beispiel 1:                    Kp(s+0.2)
%                Go =  -------------------------------------
%                     (s^2 + 2s + 10)(s+4)(s^2 + 0.2s +0.1)
%   
figure(8);
z=[1 6.2 19.7 44.2 109.8 24];
n=[1 6.2 19.7 44.2 9.8 4];
System8=tf(z, n);
O=logspace(-2, 1, 600);
[Real, Imag]=nyquist(System8, O);
plot(Real(:), Imag(:), 'b', Real(:), -Imag(:), 'b--');
hold on;
plot([0 12 11.6],[0 -10 -9], 'b', [12 11], [-10 -9.7], 'b')
axis([-5 21 -13 13]);
xticks([-5 0 5 10 15 20]);
yticks([-10 -5 0 5 10]);
axis('square');
grid('on');
latextitle('Ortskurve der R\"uckf\"uhrdifferenzfunktion $$(k_P=100)$$');
plot(Real(1), 0, 'bo');
plot(1, 0, 'bx');
latextext(5.5, 2.5, '$$\omega=0$$');
latextext(0.5, 2.5, '$$1$$');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextext(7, -5, '$$F(s)$$');
latextext(7.5, -12, '$$\leftarrow \omega$$');
hold off;


%%  Beispiel 3 mit Totzeit  :           1
%                             Go = ---------- exp(-sTt)
%                                  s(s+1)(s+3)
%   
figure(9);
z=[1];
n=[1 4 3 0];
System9=tf(z, n);
W=logspace(-1, 0.1, 300);
[Betrag, Phase]=bode(System9, W);    % Phase in Grad
k=0;
for Tt=[0, 1, 2, 3, 4]
    PhaseT=Phase(:)*pi/180-Tt*W';
    Real=Betrag(:).*cos(PhaseT);
    Imag=Betrag(:).*sin(PhaseT);
    plot(Real(:), Imag(:), 'b');
    k=k+1;
    st(k,1)=Real(1,1);
    hold on
end
plot([-10 10], [0 0], 'k:', [0 0], [-20 10], 'k:');
plot([-1 -1],[0.1 -0.1],'k', [-1.1 -0.9],[0 0],'k');
%   Kreis
omega=0.18;                     %  Berechnung des Ortskurvenpunktes für diese Frequenz
[B, P]=bode(System9, omega);
plot(B*cos(P*pi/180), B*sin(P*pi/180), 'bx');
plot(B*cos(2*pi/360*[P:-1:-160]), B*sin(2*pi/360*[P:-1:-160]), 'b--');
%grid('on');
latextitle('Offene Kette mit I-Verhalten und Totzeit');
axis([-3 1.6 -3.8 0.8]);
xticks([-3 -2 -1 0 1]);
yticks([-3 -2 -1 0 1]);
axis('square');
latextext(-2.6, -3.55, '$$T_{\rm t}=$$');
latextext(st(1,1)-0.1, -3.5, '0');
latextext(st(2,1)-0.1, -3.5, '1');
latextext(st(3,1)-0.1, -3.5, '2');
latextext(st(4,1)-0.1, -3.5, '3');
latextext(st(5,1)-0.1, -3.5, '4');
latextext(-2.5, -2.5, '$$\omega \uparrow$$')
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold off




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