%% ...
%% ... Example from Schmidt: 
%% ... The engineering of chemical reactions, Oxford University Press, Inc
%% ... 1998 
%% ...
%% ... Milan Bernauer, ICT, Prague 2013
%% ...
clc; clear all;
tb = 20;            % [min]
xa = 0.8;           % [1]
k  = -log(1-xa)/tb; % [1/min]
n  = 100;           %   
ti = zeros(1,n);
tt = zeros(1,n);
xi = zeros(1,n);
xt = zeros(1,n); 
tii= (1./(1-xa).^(1./n) - 1) ./ k;

for i = 1: n
    ti(i) = (1./(1-xa).^(1./i) - 1) ./ k;
    tt(i) = i.*ti(i);
    xt(i) = 1 - (1+k.*tii).^(-i);
    xi(i) = (1+k*tii).^(-(i-1)).*(1-1./(1+k*tii));
end

subplot(2,1,1)
yrange=[0 max(ti)];
xrange=[0 n];
limits=[xrange yrange];
plot((1:n),ti,'-',(1:n),tt,'--')
title(['Residence time for n = 1 to ',num2str(n),' CSTRs for X_A =',num2str(xa*100) ,'%'])
grid on
axis(limits)
xlabel('N of CSTR')
ylabel('t_R  / min')
legend('t_R(i-th CSTR)','t_R(Total)')

subplot(2,1,2)
yrange=[0 max(xi)];
xrange=[0 n];
limits=[xrange yrange];
[AX,H1,H2] = plotyy((1:n),xi,(1:n),xt);
title(['Total conversion and conversion in every (single) element of serie of ',num2str(n),' CSTRs'])
grid on
axis(limits)
xlabel('CSTR(i)')
set(get(AX(1),'Ylabel'),'String','X_A(i)') 
set(get(AX(2),'Ylabel'),'String','X_A(1 to i)') 
