function PFRlength
% 等溫管式反應器的管長計算
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/02/21 $

clear all
clc

global PA0 eps

d = 0.1;                % 反應器內經，m
F0 = 10;                % 進氣速率，kmol/h
P = 1;                  % 反應壓力，atm
yA0 = 2/3;              % A的初始摩爾分率
PA0 = yA0 * P;          % A的初始分壓，atm
delta = -1/2;
eps = yA0 * delta;      % eps: 膨脹因子
xA = [0.1:0.1:0.5];     % A的摩爾轉化率

% 計算不同轉化率下的管長
for i = 1:length(xA)
    I = quadl(@func,0,xA(i));   % 用自適應Lobatto法求積,積分限為[0,xA]
    L(i) = 4*yA0*F0/(pi*d^2)*I; % L: 管長,m
end
xAe = fzero(@Rate,0.5);         % 平衡轉化率

% 結果顯示
fprintf('Results:\n\txA:\t')
fprintf('\t%.2f',xA)
fprintf('\n\tL (m):')
fprintf('\t%.2f',L)
fprintf('\n\txAe = %.3f',xAe)

% ------------------------------------------------------------------
function f = func(xA)  
f = 1./Rate(xA);

% ------------------------------------------------------------------
function r = Rate(xA)
global PA0 eps
PA = PA0 * (1-xA)./(1+eps*xA);
PS = PA0/2 * xA./(1+eps*xA);
r = 116 * (PA.*PA - PS/1.27);