%   Script KOPPELM.M zur Herstellung von Bildern
%
%   Berechnung der Koppelmaße für die Klärschlammverbrennungsanlage
%
% verwendet EPSBILD.M
%
% J. Lunze
% 14.10.1996
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 25.3.2004
% für 5. Auflage: 18.10.2007
% für 7. Auflage 31.3.2012
% für 8. Auflage 3.4.2014
% für 9. Auflage 9.5.2016
% für 10. Auflage: 20.10.2019
%
%
clear
close all
format short e
%  statische Verstärkungsmatrix für ein System mit 2 Ausgängen und 5 Eingängen:
r=2;
m=5;
Ks= [0.2  454.8  -31.6  8.8  2.543
     0    -44.88  -2.07 2.12  0];
%%
%
%   Koppelmaß KKs aus Statikmatrix
%
%
for i=1:2
for j=1:5
    KKs(i,j) = Ks(i,j)*Ks(i,j);
end;
end;
kmax=max(max(KKs));
KKs=KKs/kmax

%%
%
%   Koppelmaß KKd aus Gewichtsfunktion
%
%
%   Beschreibung jedes Elementes durch ein PT1-Glied
%   Zeitkonstanten in T1:
T1=[0  1125  0  902.8  0
    0  131   0  127.8  0];
KKd=zeros(r,m);
for i=1:2
for j=[2, 4]
Tend=T1(i,j)*10;
ke=1000;
dt=Tend/ke;
I=0;
for k=1:ke
    I=I+exp(-2*k*dt/T1(i,j));
end;
I=I*dt/T1(i,j)/T1(i,j);
KKd(i,j)=I;
end;
end;
kmax=max(max(KKd));
KKd=KKd/kmax

%%   Gewichtsfunktionsmatrix der Klärschlammverbrennungsanlage
%
%
figure(1);
T=[0:10:4000];
Gewicht=tf([Ks(1,4)],[T1(1,4) 1]);
[y12, T, x12]=impulse(Gewicht, T);
Gewicht2=tf([Ks(1,2)],[T1(1,2) 1]);
[y11, T, x11]=impulse(Gewicht2, T);
Gewicht3=tf([Ks(2,2)],[T1(2,2) 1]);
[y21, T, x21]=impulse(Gewicht3, T);
Gewicht4=tf([Ks(2,4)],[T1(2,4) 1]);
[y22, T, x22]=impulse(Gewicht4, T);
subplot(2,2,1);
plot(T, y11, 'b', [0 4000], [0 0], 'k:');
xleer;
latexylabel('$$g_{11}$$');
rotateY;
yticks([0 0.4]);
latextitle('Gewichtsfunktionsmatrix');
axis([0 4000 -0.1 0.5]);
subplot(2,2,2);
plot(T, y12, 'b', [0 4000], [0 0], 'k:');
xleer;
latexylabel('$$g_{12}$$');
rotateY;
yticks([0 0.01]);
latextitle('Kl\"arschlammverbrennung');
axis([0 4000 -0.002 0.01]);
subplot(2,2,3);
plot(T, y21, 'b', [0 4000], [0 0], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$g_{21}$$');
rotateY;
yticks([-0.4 0]);
axis([0 4000 -0.4 0.05]);
xticks([0 2000 4000]);
%
subplot(2,2,4);
plot(T, y22, 'b', [0 4000], [0 0], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$g_{22}$$');
rotateY;
yticks([0 0.02]);
xticks([0 2000 4000]);
axis([0 4000 -0.005 0.02]);


%%   normierte Gewichtsfunktionsmatrix
%     der Klärschlammverbrennungsanlage
%
%
figure(2);
T=[0:10:4000];
Klaerschlamm=tf([1/2/T1(1,4)],[T1(1,4)/2 1]);
[y12, T, x12]=impulse(Klaerschlamm, T);
Klaerschlamm2=tf([1/2/T1(1,2)],[T1(1,2)/2 1]);
[y11, T, x11]=impulse(Klaerschlamm2, T);
Klaerschlamm3=tf([1/2/T1(2,2)],[T1(2,2)/2 1]);
[y21, T, x21]=impulse(Klaerschlamm3, T);
Klaerschlamm4=tf([1/2/T1(2,4)],[T1(2,4)/2 1]);
[y22, T, x22]=impulse(Klaerschlamm4, T);
subplot(2,2,1);
plot(T, y11, 'b', [0 4000], [0 0], 'k:');
xleer;
latexylabel('$$\left| \frac{g_{11}}{k_{\rm s11}}\right|^2$$');
rotateY;
yticks([0 0.00006]);
latextitle('Normierte ');
axis([0 500 0 0.00006]);
subplot(2,2,2);
plot(T, y12, 'b', [0 4000], [0 0], 'k:');
xleer;
latexylabel('$$\left| \frac{g_{12}}{k_{\rm s12}}\right|^2$$');
rotateY;
yticks([0 0.00006]);
latextitle('Gewichtsfunktionsmatrix');
axis([0 500  0 0.00006]);
subplot(2,2,3);
plot(T, y21, 'b', [0 4000], [0 0], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$\left| \frac{g_{21}}{k_{\rm s21}}\right|^2$$');
rotateY;
yticks([0 0.00006]);
axis([0 500 0 0.00006]);
xticks([0 200 400]);
%
subplot(2,2,4);
plot(T, y22, 'b', [0 4000], [0 0], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$\left| \frac{g_{22}}{k_{\rm s22}}\right|^2$$');
rotateY;
yticks([0 0.00006]);
xticks([0 200 400]);
axis([0 500 0 0.00006]);

%%  Kopplungsanalyse mit Koppelfaktor kappa
%
%
figure(3);
subplot(4,1,[1,2]);
W=logspace(-5, 1, 100);
Koppelfaktor=tf([Ks(1,4)],[T1(1,4) 1]);
[G12, P12]=bode(Koppelfaktor,W);
Koppelfaktor2=tf([Ks(1,2)],[T1(1,2) 1]);
[G11, P11]=bode(Koppelfaktor2,W);
Koppelfaktor3=tf([Ks(2,2)],[T1(2,2) 1]);
[G21, P21]=bode(Koppelfaktor3,W);
Koppelfaktor4=tf([Ks(2,4)],[T1(2,4) 1]);
[G22, P22]=bode(Koppelfaktor4,W);
Gsize=max(size(G12));
for k=1:Gsize
   Kappa(k)=G12(k)/G11(k)*G21(k)/G22(k);
   Phi(k)=P12(k)+P21(k)-P11(k)-P22(k);
   Kkompl(k)=Kappa(k)*cos(Phi(k)*pi/180) + Kappa(k)*sin(Phi(k)*pi/180)*i;
end;
semilogx(W, Kappa, 'b');
latextitle('Koppelfaktor f\"ur die Kl\"arschlammverbrennungsanlage');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
latexylabel('$$|\kappa|$$');
rotateY;
yticks([0.2 0.6]);
xticks([0.00001 0.001 0.1]);
axis([0.00001  1 0.2 0.6]);

%%  Relative Gain Array
%
%
figure(4);
for k=1:Gsize
  RGA11(k)=abs(1/(1-Kkompl(k)));
  RGA12(k)=Kappa(k)*RGA11(k);
end;
subplot(2,2,1);
semilogx(W, RGA11, 'b');
hold on;
plot([0.1 0.1], [0 1], 'k:', [0.5 0.5], [0 1], 'k:');
axis([0.00001, 10, 0, 1]);
xleer;
latexylabel('$$RGA_{11}$$');
rotateY;
yticks([0 1]);
latextitle('RGA ');
hold off;
subplot(2,2,2);
semilogx(W, RGA12, 'b');
hold on;
plot([0.1 0.1], [0 1], 'k:', [0.5 0.5], [0 1], 'k:');
xleer;
latexylabel('$$RGA_{12}$$');
rotateY;
yticks([0 1]);
latextitle('Kl\"arschlammverbrennung')
axis([0.00001, 10, 0, 1]);
hold off;
subplot(2,2,3);
semilogx(W, RGA12, 'b');
hold on;
plot([0.1 0.1], [0 1], 'k:', [0.5 0.5], [0 1], 'k:');
latexxlabel('$$\omega$$');
latexylabel('$$RGA_{21}$$');
rotateY;
yticks([0 1]);
axis([0.00001, 10, 0, 1]);
hold off;
subplot(2,2,4);
semilogx(W, RGA11, 'b');
hold on;
plot([0.1 0.1], [0 1], 'k:', [0.5 0.5], [0 1], 'k:');
latexxlabel('$$\omega$$');
latexylabel('$$RGA_{22}$$');
rotateY;
yticks([0 1]);
hold off;
axis([0.00001, 10, 0, 1]);

%%   Ausgabe der Bilder als .eps-Dateien
%
Epsbildc15('koppelm');
