function Absorption
% 氣體吸收模擬計算
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/20 $

clear all
clc

global alpha C0 C1 y7 L H
H = 5;      % liquid holdup，kg
L = 100;    % kg/h
C0 = 0.025;
C1 = 0.05;
y7 = 0;     % 不是0.002;
G = 500;    % kg/h
k = 0.4;
alpha = k*G/L;
dt = 1/120;
stoptime = 11*dt;
[t,y] = ode45(@MassBalanceEqs,[0:dt:stoptime],[0 0 0 0 0 0]);
tr = spline(y(:,6),t,0.002);  % tr: the time required, h

% Output results
disp('Results (t — hour, yi — wt%):')
disp('       t        y1        y2        y3        y4        y5        y6')
disp([t y*100])
fprintf('\tTime required: %.1f (min)', tr*60)

% ------------------------------------------------------------------
function dydt = MassBalanceEqs(t,y)
% ODEs -- Mass Balances for i-th plate of absorption column (i=1~6)
global alpha C0 C1 y7 L H
if t==0
    y0 = C0;
end
if t>0 
    y0 = C1;
end
% i=1
dydt(1) = (alpha*y0-(alpha+1)*y(1)+y(2))*L/H;
% i=2...5
for i=2:5
    dydt(i) = (alpha*y(i-1)-(alpha+1)*y(i)+y(i+1))*L/H;
end
% i=6
dydt(6) = (alpha*y(5)-(alpha+1)*y(6)+y7)*L/H;
dydt = dydt';