function PFRconv
% 已知等溫管式反應器的管長,計算出口轉化率
%
%   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
L = 5;                  % 反應管長，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 = fzero(@Func4Fzero,0.3,[],d,L,F0,yA0);         % 轉化率
fprintf('所求的轉化率為：\txA = %.4f',xA)

% ------------------------------------------------------------------
function f = Func4Fzero(xA,d,L,F0,yA0)  
I = quadl(@Func,0,xA);      % 用自適應Lobatto法求積,積分限為[0,xA]
f = L - 4*yA0*F0/(pi*d^2)*I;

% ------------------------------------------------------------------
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);