%   Script BOBREAK.M zur Herstellung von Bildern
%
%   Beobachtbarkeit zweier R\"uhrkesselreaktoren
%
%  verwendet BSPREA3.M, EPSBILD.M
%
% J. Lunze
% 3.10.1996
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 25.3.2004
% für 5. Auflage: 13.1.2008
% für 7. Auflage: 7.5.2012
% für 10. Auflage: 3.9.2019
%
%
clear
close all
Bsprea3;     % definiert Reaktor=ss(A, B, C, D);
%%
%
%     E/A-Verhalten der R\"uhrkesselreaktoren
%
%
figure(1);
Te=3;
dt=0.005;
T=[0:dt:Te];
U=2*[0:dt:Te];
kmax=Te/dt+1;
x0=[2; 4];
[Y, T, X]=lsim(Reaktor, U, T, x0);
subplot(6,1,[1,2]);
plot(T, U,'b');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latextitle('Verhalten zweier R\"uhrkesselreaktoren');
xleer;
axis([0 Te 0 2*Te]);
yticks([0 3 6]);
subplot(6,1,[3,4]);
plot(T, Y,'b');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latexxlabel('$$t$$ in min');
yticks([0 4]);
axis([0 Te 0 4]);
xticks([0 1 2 3 4]);
yticks([0 2 4]);



%%    Berechnung der freien Bewegung
%
%
figure(2);
subplot(4,1,[1,2]);
[Yerzw, T,  X]=lsim(Reaktor, U, T, [0;0]);
Y=Y-Yerzw;
plot(T, Y,'b');
hold on;
latexylabel('$$y_{\rm frei}$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
yticks([1.26 2.62]) ;
latexxlabel('$$t$$ in min');
xticks([0.5 2]);
latextitle('Eigenbewegung der R\"uhrkesselreaktoren');
axis([0 Te 0 4]);
%     Zwei Meßpunkte herausgreifen bei t1=0,5 und t2=2
t1=0.5;
t2=2;
k1=t1/dt+1;
plot(t1, Y(k1,1), 'bo');
k2=t2/dt+1;
plot(t2, Y(k2,1), 'bo');
plot([t1 t1], [0 Y(k1,1)], 'k:', [0 t1], [Y(k1,1) Y(k1,1)], 'k:');
plot([t2 t2], [0 Y(k2,1)], 'k:', [0 t2], [Y(k2,1) Y(k2,1)], 'k:');
hold off;


%%   Berechnung des Anfangszustandes
%
%
%
%        Berechnung von W(Te)
%
dt1=Te/1000;
CC=C'*C;
expA=expm(A*Te);
W=CC+expA'*CC*expA;
W=W/2;
for k=1:Te/dt1-1
  t=k*dt1;
  expAt=expm(A*t);
  W=W + expAt'*CC*expAt;
end;
W=W*dt1;
%
Int = C'*Y(1,1) + expA'*C'*Y(kmax,1);
for k=1:kmax
  t=k*dt;
  expAt=expm(A*t);
  Int=Int + expAt'*C'*Y(k,1);
end;
Int=Int*dt;
x0=inv(W)*Int




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