function Semi_BatchReactor
% 半連續反應器—苯的氯化
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/16 $

clear all
clc
global k1 k2 k3 V nB0 F
k1 = 510;       % [m3/kmol.h]
k2 = 64;        % [m3/kmol.h]
k3 = 2.1;       % [m3/kmol.h]
nB0 = 50;       % kmol
V = 1.46*50;    % m3
F = 1.44*nB0;
n0 = [nB0 0 0 0 0];
[t,y] = ode45(@MassBalances,[0:0.01:2.0],n0);
tOpt4M = t(find(y(:,2)==max(y(:,2))))
tOpt4D = t(find(y(:,3)==max(y(:,3))))
plot(t,y(:,1),'k-',t,y(:,2),'b--',t,y(:,3),'k-.',t,y(:,4),'r--',t,y(:,5),'b-')
xlabel('Time (hr)')
ylabel('Concentration (moles)')
legend('n_B','n_M','n_D','n_T','n_C')
str = {'B','M','D','T','C'};
for i=1:5
    gtext(str{i})
end
fprintf('\t氯苯達到最大產量所對應的最優時間為: %.2f %s\n',tOpt4M,'hr')
fprintf('\t二氯苯達到最大產量所對應的最優時間為: %.2f %s',tOpt4D,'hr')

% ------------------------------------------------------------------
function dydt = MassBalances(t,y)
global k1 k2 k3 V nB0 F
nB = y(1);
nM = y(2);
nD = y(3);
nT = y(4);
nC = y(5);
if nC > 0.12*nB0
    nC = 0.12*nB0;
    dnCdt = 0;
else
    dnCdt = F - k1*nB*nC/V - k2*nM*nC/V - k3*nD*nC/V;
end
dnBdt = -k1*nB*nC/V;
dnMdt = k1*nB*nC/V - k2*nM*nC/V;
dnDdt = k2*nM*nC/V - k3*nD*nC/V;
dnTdt = k3*nD*nC/V;
dydt = [dnBdt; dnMdt; dnDdt; dnTdt; dnCdt];