function AbsPackedHeight
% Gas Absorption in Packed Towers
% 計算填料吸收塔的總傳質單元數NOG
%
% 采用自適應Lobatto求積(Adaptive Lobatto quadrature)函數quadl()進行數值積分
% 
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/02/23 $

clear all
clc

% 離散平衡數據
xi = 1e-3*[1.96  1.40	0.846	0.562	0.422	0.281   0.141   0.056];
yi = 1e-3*[51.3	 34.2	18.6	11.2	7.63	4.21    1.58    0.658];

G2L = 1/40;     % 氣液比   
x2 = 0;         % 塔頂液相摩爾分率    
y2 = 0.005;     % 塔頂氣相摩爾分率
y1 = 0.055;     % 塔底氣相摩爾分率
HOG = 0.69;     % 傳質單元高度(m)

sp = csaps(xi,yi,1);    % 與sp = spline(xi,yi)等效

% 畫擬合曲線，直觀地檢查擬合效果是否良好
fnplt(sp); hold on,plot(xi,yi,'ok'), hold off

NOG = quadl(@func,y2,y1,[],[],sp,G2L,x2,y2);  % 總傳質單元數NOG

disp('計算結果:')
fprintf('\n   總傳質單元數為：%.2f\n',NOG)
fprintf('\n   填料塔高度為：%.2f %s\n',NOG * HOG,'（m）')

% ------------------------------------------------------------------
function f = func(y,sp,G2L,x2,y2)    % 定義被積函數
var = x2/(1-x2) + G2L*(y./(1-y)-y2/(1-y2));
x = var./(1+var);           % 操作方程
yeq = fnval(sp,x);          % 平衡關系
f = 1./(y - yeq);