%   Script BODEDIA.M zur Erzeugung von Bildern
%
%   Bodediagramme (Geradenapproximationen)
%
% J. Lunze
% 23.2.1996
% für 2. Auflage: 14.12.1998
% für 3. Auflage: 3.4.2001
% für 4. Auflage: 10.9.2003
% für 5. Auflage: 18.4.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 28.4.2008
% für 8. Auflage: 5.1.2010
% für 9. Auflage: 13.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 13.9.2019
%
%   
clear

close all
%set(0,'DefaulttextFontname','times');
%set(0,'DefaulttextFontsize',24);

%%   Geradenapproximation des Amplitudenganges
% 
figure(1);
latextitle('Geradenapproximationen des Amplitudenganges');
W=logspace(-1, 1, 200);
subplot(3,2,1);
z=[1];
n=[1 0];
System=tf(z, n);
[B, P, W]=bode(System, W);
semilogx(W, 20*log10(B(:)), 'b');
xleer;
latextitle('Amplitudeng\"ange');
latexylabel('$$|G|$$ in dB');
rotateY;
axis([0.1 10 -20 20]);
yticks([-20 0 20]);
%
subplot(3,2,2);
z=[1 0];
n=[0.001 1];
System=tf(z, n);
[B, P, W]=bode(System, W);
semilogx(W, 20*log10(B(:)), 'b');
xleer;
axis([0.1 10 -20 20]);
yticks([-20 0 20]);
%
subplot(3,2,3);
z=[1];
n=[1 1];
System=tf(z,n);
[B, P, W]=bode(System, W);
semilogx(W, 20*log10(B(:)), 'b');
xleer;
hold on;
plot([0.1 1 10],[0 0 -20],'b--');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;
axis([0.1 10 -20 20]);
yticks([-20 0 20]);
%
subplot(3,2,4);
z=[1 1];
n=[0.001 1];
System=tf(z, n);
[B, P, W]=bode(System, W);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
plot([0.1 1 10],[0 0 20],'b--');
xleer;
hold off;
axis([0.1 10 -20 20]);
yticks([-20 0 20]);
%
subplot(3,2,5);
z=[1];
n=[1 0.4 1];
System=tf(z,n);
[B, P, W]=bode(System, W);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
plot([0.1 1 10],[0 0 -40],'b--');
latexylabel('$$|G|$$ in dB');
rotateY;
latexxlabel('$$\omega$$ in rad/s');
hold off;
axis([0.1 10 -20 20]);
yticks([-20 0 20]);
xticks([0.1 1 10]);
%
subplot(3,2,6);
z=[1 0.4 1];
n=[0.0001 0.01 1];
System=tf(z, n);
[B, P, W]=bode(System, W);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
plot([0.1 1 10],[0 0 40],'b--');
latexxlabel('$$\omega$$ in rad/s');
hold off;
axis([0.1 10 -20 20]);
yticks([-20 0 20]);
xticks([0.1 1 10]);
 

%%  |G(jwo)| und Resonanzüberhöhung in Abhängigkeit von
%   der Dämpfung
%   
figure(2);
k=0;
for d=0.01:0.02:1/sqrt(2)
   k=k+1;
   D1(k,1)=d;
   D2(k,1)=d;
   Go(k,1)=20*log10(1/2/d);
   Gr(k,1)=20*log10(1/(2*d*sqrt(1-d*d)));
end
for d=1/sqrt(2):0.02:0.99
   k=k+1;
   D2(k,1)=d;
   Go(k,1)=20*log10(1/2/d);
end
plot(D2, Go, 'b', D1, Gr, 'b--');
grid on;
latextitle('Resonanz\"uberh\"ohung');
latexxlabel('$$d$$');
latexylabel('\begin{minipage}{1.5cm}$$|G_0|, |G_{\rm r}|$$ in dB\end{minipage}');
rotateY;
axis([0 1 -10 30]);
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([-10 0 10 20 30]);
 

%%     Schwingfähige Systeme
%
%   
figure(3);
T=0.5;
d=0.1;
zs=1;
ns=[T*T 2*d*T 1 0];
System=tf(zs, ns);
pole=max(abs(roots(ns)));
W=logspace(-1, 1, 1000);
[B, P]=bode(System, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
Bwmin=20*log10(B(1));
Wmin=W(1,1);
Wmax=max(W);
Bpz=Bwmin-20*log10(pole/Wmin);
Bwmax=Bpz-60*log10(Wmax/pole);
axis([0.1 10 -40 20]);
yticks([-40 -20 0 20]);
plot([Wmin pole Wmax], [Bwmin Bpz  Bwmax],'b--');
plot([Wmin Wmax], [0 0], 'k:');
plot([pole pole], [20 -40], 'k:');
latextitle('Bodediagramm eines schwingf\"ahigen Systems');
xleer;
%xticks([0.1 1 10]);
grid('on');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;
%
subplot(2, 1, 2);
semilogx(W, P(:), 'b');
hold on
plot([pole pole], [20 -300], 'k:');
axis([0.1 10 -300 0]);
latexxlabel('$$\omega$$ in rad/s');
xticks([0.1 1 10]);
latexylabel('$$\phi$$ in Grad');
rotateY;
yticks([-270 -180 -90 0]);
grid('on');
 
%%
figure(4);
zs=1.6*[0.7 0 1];
ns=[0.1 0 1 0 0];
System=tf(zs, ns);
pole=max(abs(roots(ns)));
zero=max(abs(roots(zs)));
W=logspace(-1, 1, 1000);
[B, P]=bode(System, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(B(:)), 'b');
grid on;
hold on;
Bwmin=20*log10(B(1));
Wmin=W(1,1);
Wmax=max(W);
Bpz=Bwmin-40*log10(zero/Wmin);
Bwmax=Bpz-40*log10(Wmax/pole);
axis([0.1 10 -40 40]);
yticks([-40 -20 0 20 40]);
plot([Wmin zero pole Wmax], [Bwmin Bpz Bpz Bwmax],'b--');
latextext(0.9, 30, 'Plateau');
latextext(1.5, 10, '$$\downarrow$$');
latextitle('Bodediagramm eines schwingf\"ahigen Systems');
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;
%
subplot(2, 1, 2);
semilogx(W, P(:), 'b');             % bei MATLAB 2007:  -P(:) drucken
grid on;
axis([0.1 10 -200 20]);
latexxlabel('$$\omega$$ in rad/s');
xticks([0.1 1 10]);
latexylabel('$$\phi$$ in Grad');
rotateY;
yticks([-270 -180 -90 0]);
 
%%
figure(5);
zs=1.6*[0.41 0 1];
ns=[0.33 0 1 0 0];
System=tf(zs, ns);
pole=max(abs(roots(ns)));
zero=max(abs(roots(zs)));
W=logspace(-1, 1, 1000);
[B, P]=bode(System, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(B(:)), 'b');
grid on;
hold on;
Bwmin=20*log10(B(1));
Wmin=W(1,1);
Wmax=max(W);
Bpz=Bwmin-40*log10(zero/Wmin);
Bwmax=Bpz-40*log10(Wmax/pole);
axis([0.1 10 -40 40]);
yticks([-40 -20 0 20 40]);
plot([Wmin zero pole Wmax], [Bwmin Bpz Bpz Bwmax],'b--');
latextext(0.5, 30, 'Plateau');
latextext(1, 10, '$$\searrow$$');
latextitle('Bodediagramm eines schwingf\"ahigen Systems');
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;
%
subplot(2, 1, 2);
semilogx(W, P(:), 'b');             % bei MATLAB 2007:  -P(:) drucken
grid on;
axis([0.1 10 -200 20]);
latexxlabel('$$\omega$$ in rad/s');
xticks([0.1 1 10]);
latexylabel('$$\phi$$ in Grad');
rotateY;
yticks([-270 -180 -90 0]);
 

%%    Bodediagramm und Ortskurve von PTn-Gliedern
%
figure(6)
PT1=tf([1],poly(-1));
PT2=tf(1, poly([-1, -1]));
PT3=tf(1, poly([-1, -1, -1]));
PT4=tf(1, poly([-1, -1, -1, -1]));
[R1, I1]=nyquist(PT1);
[R2, I2]=nyquist(PT2);
[R3, I3]=nyquist(PT3);
[R4, I4]=nyquist(PT4);
subplot(2, 4, [1 2 5 6]);
plot(R1(:), I1(:), 'b', R2(:), I2(:), 'b', R3(:), I3(:), 'b', R4(:), I4(:), 'b');
hold on
plot([-2 2],[0 0],'k:', [0 0],[-2 2],'k:', 1, 0, 'bo');
axis([-0.4 1.1 -1 0.5]);
axis('square');
xticks([0 1]);
yticks([-1 0]);
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Ortskurven von PT$$_n$$-Gliedern');
latextext(0.5, -0.4, '$$\leftarrow$$');
hold off
%
subplot(2, 4, [3 4 7 8]);
W=logspace(0, 5, 500);
[R1, I1]=nyquist(PT1, W);
[R2, I2]=nyquist(PT2, W);
[R3, I3]=nyquist(PT3, W);
[R4, I4]=nyquist(PT4, W);
plot(R1(:), I1(:), 'b', R2(:), I2(:), 'b', R3(:), I3(:), 'b', R4(:), I4(:), 'b');
hold on
plot([-2 2],[0 0],'k:', [0 0],[-2 2],'k:');
axis([-0.0003 0.0003 -0.0003 0.0003]);
axis('square');
xticks([-0.0002 0.0002]);
latexxlabel('Re');
%latexylabel('Im');
yleer
latextext(0.1e-4, -2e-4, '$$\uparrow  {\rm PT}_1$$');
latextext(-2.4e-4, -0.5e-4, '$${\rm PT}_2  \rightarrow$$');
latextext(-2e-4, 1e-4, '$${\rm PT}_3  \downarrow$$');
latextext(1e-4, 0.2e-4, '$$\leftarrow {\rm PT}_4$$');
hold off



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