function DistStagesCal
% DistStagesCal.m
% 計算雙組分簡單精餾塔的理論板數
%
% 主程序采用MATLAB的自適應Lobatto求積(Adaptive Lobatto quadrature)函數quadl()
% 進行數值積分；
% 被積函數func1()和func2()中利用三次樣條擬合法由已知離散數據確定y與x的函數關系
% 
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/02/22 $

clear all
clc

xi = [0.178	0.275	0.372	0.456	0.650	0.844]; % 液相平衡濃度
yi = [0.243	0.382	0.518	0.616	0.795	0.931]; % 汽相平衡濃度
xf = 0.4;   % 進料組成
xd = 0.9;   % 塔頂組成
xw = 0.15;  % 塔底組成
R = 5;      % 精餾段回流比   
R1 = 4;     % 提餾段回流比

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

% 畫擬合曲線，直觀地檢查擬合效果是否良好
x4plot = linspace(xi(1),xi(end),200);
y4plot = fnval(sp,x4plot);
plot(xi,yi,'o',x4plot,y4plot,'-')

N = quadl(@func1,xf,xd,[],[],sp,xd,R);  % 精餾段理論板數N
M = quadl(@func2,xw,xf,[],[],sp,xw,R1); % 提餾段理論板數M
N = round(N+0.5);       % 圓整
M = round(M+0.5);

disp('計算結果:')
fprintf('\n   精餾段理論板數N為：%d %s\n',N,'（塊）')
fprintf('   提餾段理論板數M為：%d %s\n',M,'（塊）')
fprintf('   共需理論板數為：%d %s\n',N+M,'（塊）')

% ------------------------------------------------------------------
function f = func1(x,sp,xd,R)        % 定義被積函數(精餾段)
y = fnval(sp,x);
f = 1./(y-x-(xd-y)/R);

% ------------------------------------------------------------------
function f = func2(x,sp,xw,R1)         % 定義被積函數(提餾段)
y = fnval(sp,x);
f = 1./(y-x-(y-xw)/R1);