%   Script PTTBODE.M
%
%   Bodediagramm eines Totzeitsystems
%
% J. Lunze
% 12.12.1995
% für 2. Auflage: 14.12.1998
% für 3. Auflage: 3.4.2001
% für 4. Auflage: 15.9.2003
% für 5. Auflage: 20.6.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 23.5.2008
% für 9. Auflage: 15.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 23.03.2020
%
%   
clear
close all
z=[1];


%%     Bodediagramm des Totzeitsystems exp(-sT) für T=1
%
figure(1);
Tt=1;
W=logspace(-1, 1, 100);
for k=1:length(W)
    Betrag(k,1)=0;
    Phase(k,1)=-(W(k)*Tt)*360/2/pi;   % Umrechnung in Grad
end;
subplot(2, 1, 1);
semilogx(W, Betrag, 'b');
%axis([0.1 10 -50 50]);
yticks([-1 0 1]);
latexylabel('$$|G|$$ in dB');
rotateY;
yticks([0]);
xleer;
grid('on');
latextitle('Totzeitglied $$G(s) = e^{-sT_{\rm t}}$$');
subplot(2, 1, 2);
semilogx(W, Phase, 'b');
hold on
grid('on');
latexylabel('$$\phi$$ in Grad');
rotateY;
latexxlabel('$$\omega/\omega_0$$');
axis([0.1 10 -360 0]);
xticks([0.1 1 10]);
yticks([-360  -180  0]);
latextext(1.05, -280, '$$\omega_0=\frac{1}{T_{\mathrm{t}}}$$');
latextext(1.1, -40, '$$-57^o$$');
plot(1, -57, 'bo');
plot(1, -360, 'bo');


%% Approximation des Totzeitgliedes durch ein PTn-Glied
%  
figure(2);
subplot(3,1,[1,2]);
T=0:0.01:2.5;
for n=[1 2 4 7 10 40]
    A = -eye(n, n)*n;
    Nenner=poly(A);
    System=tf([n^n], Nenner);
    [Y, T]=step(System, T);
    plot(T, Y, 'b');
    hold on;
    axis([0 2 0 1]);
end;
plot([0 1 1 2.5], [0 0 1 1], 'b');
axis([0 2.5 -0.05 1.1]);
xticks([0 1 2]);
yticks([0 0.5 1]);
latextitle('PT$$_n$$-Approximationen des Totzeitgliedes');
latexylabel('$$h$$');
rotateY;
latexxlabel('$$t$$ in s');
latextext(0.4, 0.52, '$$n=$$');
latextext(0.5, 0.45, '1');
latextext(0.5, 0.34, '2');
latextext(0.5, 0.21, '4');
latextext(0.5, 0.11, '7');
latextext(0.6, 0.05, '10');
latextext(0.8, 0.05, '40');
latextext(1.05, 0.2, 'Totzeitglied');
hold off;


%%   Übergangsfunktion eines Allpassgliedes
%
figure(3);
subplot(3,1,[1,2]);
System1=tf([-1 1], [1 1]);
[Y, T, X]=step(System1);
System2=tf([1],[1 1]);
[Y1, X1]=step(System2 ,T);
plot(T, Y, 'b', [-1 0 0],[0 0 -1], 'b', T, Y1, 'b--', [0 5],[0 0],'k:',[0 0],[-1.2 1.2],'k:');
axis([-1 5 -1.1 1.1]);
xticks([-1 0 1 2 3 4 5 6]);
yticks([-1 0 1]);
latextitle('\"Ubergangsfunktion eines Allpassgliedes');
latexylabel('$$h_{\rm A}$$');
rotateY;
latextext(0.5, -0.6, 'Allpassglied');
latextext(0.3, 0.9, 'PT$$_1$$-Glied');
latexxlabel('$$t$$ in s');



%%   Ortskurve zweier Systeme mit Totzeit
%
figure(4);
T1=2;
T2=3;
T3=4;
A=diag([-1/T1 -1/T2 -1/T3]);
Nenner=poly(A);
k1=3/T1/T2/T3;
k2=1;
Tt=60;
O1 = logspace(-3, 0, 500);
O=[];
for k=1:430
  O(1,k)=O1(1,k);
end;
System3=tf(k1, Nenner);
[Betrag1, Phase1]=bode(System3,O);
Realteil1=diag(Betrag1(:))*cos(Phase1(:)*pi/180);
Imagteil1=diag(Betrag1(:))*sin(Phase1(:)*pi/180);
Betrag=k2;
phi=-O'*Tt;
Realteil2=Betrag*cos(phi);
Imagteil2=Betrag*sin(phi);
Real=Realteil1+Realteil2;
Imag=Imagteil1+Imagteil2;
subplot(2,2,2);
plot(Real, Imag, 'b', [-2 4],[0 0],'k:',[0 0],[-4 1],'k:');
latextitle('PT$$_3$$-Glied parallel PT$$_t$$-Glied');
axis([-2 3 -4 1]);
xticks([-2 0 2]);
yticks([-4 -2 0]);
axis('square');
ylabel('');
rotateY;
latexxlabel('Re');
%
W=logspace(-3, 0.5, 1500);
Tt2=10;
System4=tf(1, [1,1]);
[Betrag2, Phase2]=bode(System4,W);
phi=-W'*Tt2/pi*180;
Phase2(:)=Phase2(:)+phi(:);
Real=[];
Imag=[];
Real=diag(Betrag2(:))*cos(Phase2(:)*pi/180);
Imag=diag(Betrag2(:))*sin(Phase2(:)*pi/180);
subplot(2,2,1);
plot(Real, Imag, 'b', [-2 4],[0 0],'k:',[0 0],[-3.5 2.5],'k:');
latextitle('PT$$_1$$T$$_t$$-Glied');
axis([-1.2 1.2 -1.2 1.2]);
xticks([-1 0 1]);
yticks([-1 0 1]);
axis('square');
latexylabel('Im');
rotateY;
latexxlabel('Re');



%%   Übergangsfunktion dieser Systeme
%
T=[0:0.1:75];
System5=tf(k1, Nenner);
subplot(2,2,4)
[Y, X]=step(System5, T);
for k=1:length(T)
    if T(k) > Tt,
       Y(k,1)=Y(k,1)+k2;
    end;
end;
plot(T, Y, 'b');
axis([0 72 -0.5 4.5]);
xticks([0 20 40 60]);
yticks([0 2 4]);
ylabel('');
latexxlabel('$$t$$');
%
subplot(2,2,3);
T=0:0.05:20;
Y=step(System4, T);
k1=0;
for k=1:length(T)
   if T(k) <= Tt2
      Y1(k,1)=0;
      k1=k1+1;
   else
      Y1(k,1)=Y(k-k1,1);
   end;
end;
plot(T, Y1, 'b');
axis([0 20 -0.2 1.2]);
xticks([0 5 10 15 20]);
yticks([0 1]);
latexylabel('$$h$$');
rotateY;
latexxlabel('$$t$$');



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