% absorp.m
% absorption example - Process Dynamics
% orginal version - 29 June 93
% (c) B. Wayne Bequette
% revised - 7 June 1997
% * * * * * * * * * * * * * * * * * * * * * * * * * *
% WARNING - 3 April 1998 - If you are using version
%   5.x of MATLAB, then you need the CONTROL TOOLBOX
%   to use the step function
% * * * * * * * * * * * * * * * * * * * * * * * * * *
% aeq = gas/liquid equil coeff
% ns = number of stages
% liq = liquid molar flowrate
% vap = vapor molar flowrate
% mhold = liquid molar holdup per stage
% xfeed = liquid feed composition
% yfeed = vapor feed composition
%
% set parameters
%
   ns = 5;
   liq = 80/60;
   vap = 100/60;
   mhold = 20/3;
   xfeed = 0;
   yfeed = 0.1;
   aeq = 0.5;
%
   amat = zeros(ns,ns);
   bmat = zeros(ns,2);
%
%  calculate ratios
%
   lvratio = -(liq+vap*aeq)/mhold;
   lratio  = liq/mhold;
   vratio  = vap*aeq/mhold;
%
   amat(1,1) = lvratio;
   amat(1,2) = vratio;
%
   for i = 2:ns-1;
     amat(i,i-1) = lratio;
     amat(i,i)   = lvratio;
     amat(i,i+1) = vratio;
   end
%
   amat(ns,ns-1) = lratio;
   amat(ns,ns)   = lvratio;
%
   bmat(1,1) = lratio;
   bmat(ns,2) = vap/mhold;
%
%  calculate the steady-state
%
   xs = -inv(amat)*bmat*[xfeed;yfeed];
%
%  perform step tests using step
%
%  step change the vapor feed composition at t = 0
%
   cmat = zeros(2,ns);
   cmat(1,ns) = 1;
   cmat(2,1)  = aeq;
   dmat = zeros(2,2);
%
  [y,x,t]=step(amat,bmat,cmat,dmat,2);
%
% actual step magnitude was 0.05
%
   y = y*0.05;
   x = x*0.05;
%
%  plot the liquid composition leaving the
%  absorber (deviation variables)
%
   plot(t,y(:,1),'w')
   xlabel('t(min)')
   ylabel('x_liq')
%
   pause
%
%  plot the vapor composition leaving the absorber
%
   plot(t,y(:,2),'w')
   xlabel('t(min)')
   ylabel('y_vap')
%
   pause
%
%  plot each stage liquid composition (deviation)
   plot(t,x,'w')
   xlabel('t(min)')
   ylabel('x_stage')
   pause
%
%  plot each stage, normalized by the final change
%
   nt = length(t);
%
   for i = 1:ns;
     xscale(:,i) = x(:,i)/x(nt,i);
   end
%
   plot(t,xscale,'w')
   xlabel('t(min)')
   ylabel('x_scaled')
   pause
%  now, do a step change in liquid feed comp
%
  [y1,x1,t1]=step(amat,bmat,cmat,dmat,1);
%
%  actual input change was 0.05
   y1 = y1*0.025;
   x1 = x1*0.025;
%
%  plot the liquid composition leaving the
%  absorber (deviation variables)
%
   plot(t1,y1(:,1),'w')
   xlabel('t(min)')
   ylabel('x_5')
%
   pause
%
%  plot the vapor composition leaving the absorber
%
   plot(t1,y1(:,2))
%
   pause
%  plot each stage liquid composition (deviation)
   plot(t1,x1,'w')
   xlabel('t(min)')
   ylabel('x_stage')
   pause
%
% plot scaled liquid compositions
%
   nt1 = length(t1);
%
   for i = 1:ns;
     xscale1(:,i) = x1(:,i)/x1(nt1,i);
   end
%
   plot(t1,xscale1,'w')
   xlabel('t(min)')
   ylabel('x_scaled')