%-----------------------------------------------------------------------
% Find the number of stages needed in binary distillation where the following
% partition function relates the vapor phase composition (y) to the liquid
% phase composition (x)
%   y(i) = alpha*x(i) / (1+(alpha-1)*x(i))
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Print out a program header -------------------------------------------
      disp('Find the number of stages needed in binary distillation')
      disp(' ')

% Input data -----------------------------------------------------------
      bottom = input('Enter the bottom composition (default=0.05): ');
      top    = input('Enter the top composition (default=0.95): ');
      alpha  = input('Enter the partition coefficient (default=1.5): ');
      if isempty(bottom) bottom = 0.05; end
      if isempty(top)    top    = 0.95; end
      if isempty(alpha)  alpha  = 1.5 ; end

% Iterate the equilibrium relationship to find the number of stages needed
      x(1) = bottom;                               %Start from bottom
      for i=1 : 99999
        y(i)   = alpha*x(i) / (1+(alpha-1)*x(i));  %Vapor phase in equilibrium with the liquid phase
        x(i+1) = y(i);                             %Vapor condenses to form liquid to the next stage
        if (x(i+1)>top)                            %Done when top composition is reached
          N = i;
          break
        end
      end

% Iterate for unknown number of times with "while" rather than "for" and "if"
% The price to pay with "while" is that we need to increment the counter index manually.
%     i = 1;
%     while x(size(x,2)) < top
%       y(i)   = alpha*x(i) / (1+(alpha-1)*x(i));  %Vapor phase in equilibrium with the liquid phase
%       x(i+1) = y(i);                             %Vapor condenses to form liquid to the next stage
%       i = i+1;
%     end
%     N = i-1;

% Print out a table ----------------------------------------------------
      disp(' ')
      disp(['The number of stages needed is ', num2str(N)])
      disp('--------------------------')
      disp(' Stage  Liquid Composition')
      disp('--------------------------')
      for i=1:N+1
        fprintf('%4i %16.5f\n', i, x(i) )
      end
      disp('--------------------------')

% Quick and dirty table ------------------------------------------------
%     [ [1:N]'  x(1:N)'  y(1:N)' ]

% Form the equilibrium line for plotting -------------------------------
      xx = 0:0.01:1;
      yy = alpha.*xx ./ (1+(alpha-1).*xx);

% Plot -----------------------------------------------------------------
% Method #1 (best). Form the steps for the McCabe-Thiele plot first, then plot
%     for i=1:N
%       u(2*i-1)=x(i); v(2*i-1)=x(i);  %Coordinates for points on the diagonal line
%       u(2*i)  =x(i); v(2*i)  =y(i);  %Coordinates for points at the corners
%     end
%     u(2*N+1)=x(N+1); v(2*N+1)=x(N+1); %One last point on the diagonal line
%     plot(xx,xx,  xx,yy,  u,v)   %Curves: diagonal, equilibrium, and steps

% Plot steps easily with "staris"; however, if we look carefully, we will
% see that there is a vertical step missing at the first stage and a
% horizontal step missing at the last stage.  For a better plot take out the
% percent signs in the last few lines and replace them with the "plot" and
% "stairs" statements

% Method #2. Let the stairs" function generate the vectors of stairs, then plot
%     [xstair ystair] = stairs(x(1:N),y(1:N));     %generate vectors of stairs
%     plot(xx,xx,  xx,yy,  xstair,ystair)

% Method #3. Plot the regular curves first, then plot the stairs on top.
      plot(xx,xx,  xx,yy),  hold on  %Diagonal and equilibrium curves
      stairs(x(1:N),y(1:N))          %Steps

% Add legends ----------------------------------------------------------
      axis([0 1 0 1])
      xlabel('Liquid Phase Composition')
      ylabel('Vapor Phase Composition')
      title('McCabe-Thiele Diagram')
      text(0.1, 0.90, ['bottom  =',num2str(bottom)])
      text(0.1, 0.85, ['top     =',num2str(top)   ])
      text(0.1, 0.80, ['# stages=',num2str(N)     ])

% Add numbers at the step corners --------------------------------------
      for i=1:2:N
        text(x(i)-0.02, y(i)+0.02, num2str(i))
      end