Code covered by the BSD License  

Highlights from
McCabeThiele

image thumbnail

McCabeThiele

by

 

10 Jun 2013 (Updated )

This function creates McCabe-Thiele Diagrams from various user inputs

mccabethiele
function mccabethiele
% Function MCCABETHIELE creates a McCabe Thiele Diagram from user inputs 
% FUNCTION MCCABETHIELE can take inputs as any of the following
% combinations: 
% External Reflux Ratio and Feed Condition
% External Reflux Ratio and External Boilup Ratio
% External Boilup Ratio and Feed Condition


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Code Written By: Vincent Cericola | ChE 212 Spring 2013 %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Updated 6/12/2013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Close Current Figures
close all

clc

% Dock Figures for Ease 
set(0,'DefaultFigureWindowStyle','docked')
f = figure;
set(f,'name','Specified Parameters','numbertitle','off')

%Change this parameter to "1" if you want the minimum steps shown in a
%separate figure window after a pause
minstep = 1;

%Change this parameter to "1" if you want the minimum reflux ratio diagram 
%made in a separate figure window after a pause
maxstep = 1;

%Change this parameter to "1" if you want the column diagram to be made in
%a separate window after a pause
column = 1;

% Print Input Choices and Corresponding Codes
fprintf('Enter the following codes to specify given parameters: \n')
fprintf('==============================================================\n')
fprintf('Enter "0" to specify reflux and feed condition \n')
fprintf('Enter "1" to specify reflux and boilup \n')
fprintf('Enter "2" to specify boilup and feed condition \n')
fprintf('==============================================================\n')


% Choose Inputs
M = input('Input Code (see above): ');
fprintf('==============================================================\n')

% Stage by Stage?
fprintf('Enter "1" to show stage steps or "0" to not show stage steps \n')
fprintf('==============================================================\n')
stages = input('Stage by Stage? :       ');
fprintf('==============================================================\n')

if stages ~= 1
    stages = 0;
    minstep = 0;
    maxstep = 0;
end

% Enter VLE Data
VLE = input('Enter VLE Data [x y] or alpha12:        ');

% Enter External Reflux Ratio
if M == 0 || M == 1
    reflux = input('Enter External Reflux Ratio:            ');
end

% Enter External Boilup Ratio
if M == 1 || M == 2
    boilup = input('Enter External Boilup Ratio:            ');
end

% Enter Feed Data
if M == 0 || M == 2
    q1 = input('Enter the value of q (feed condition):  ');
    
    if q1 ~= 1
        q = q1;
    else 
        q = 0.999;
    end
end

F = input('Enter the feed flow rate [kmol/hr]:     ');

% Enter Compositions
comps = input('Enter Xd Xb and Z [Xd Xb Z]:           ');
Xd = comps(1, 1);
Xb = comps(1, 2);
Z  = comps(1, 3);

% % To specify values instead of having user inputs, uncomment this section
% 
% % q1 = 0.7;
% % stages = 1;
% % M = 0;
% % VLE = 6;
% % VLE     = [0 0; 0.01 0.03; 0.02 0.134; 0.04 0.23; 0.06 0.304; 0.08 0.365
% %            0.1 0.418; 0.15 0.517; 0.2 0.579; 0.3 0.665; 0.4 0.729; 0.5 0.779
% %            0.6 0.825; 0.7 0.87; 0.8 0.915; 0.9 0.958; 0.95 0.979; 1 1];
% %
% % VLE = [0.0000 0.0000
% %        0.0040 0.0657
% %        0.0080 0.1380
% %        0.0160 0.2770
% %        0.0330 0.4790
% %        0.0520 0.6040
% %        0.0720 0.6750
% %        0.0940 0.7190
% %        0.1170 0.7380
% %        0.1710 0.7760
% %        0.2370 0.8000
% %        0.3180 0.8220
% %        0.4200 0.8390
% %        0.5540 0.8630
% %        0.7360 0.9090
% %        0.8000 0.9300
% %        0.9000 0.9600
% %        0.9500 0.9800
% %        1.0000 1.0000];
% % reflux  = 0.3;
% % Xd      = 0.9;
% % Xb      = 0.1;
% % Z       = 0.5;
% % q       = 0.75;
% % F       = 100;
% % boilup  = 3;

%%%% Setup VLE Curve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(VLE) > 1
    xVLE = VLE(1:size(VLE), 1);
    yVLE = VLE(1:size(VLE), 2);
    
%%%% Approximate VLE function for Stage by Stage Analysis %%%%%%%%%%%%%%%%%
    if stages == 1
        curve = 1;
        p = polyfit(xVLE,yVLE,10);
        syms K
        VLE_EQ = p(1)*K.^10 + p(2)*K.^9 + p(3)*K.^8 + p(4)*K.^7 + p(5)*K.^6 + p(6)*K.^5 + p(7)*K.^4 + p(8)*K.^3 + p(9)*K.^2 + p(10)*K+ p(11);
    end
        
else 
%%%%% Calculate VLE from Constant Relative Volatility %%%%%%%%%%%%%%%%%%%%%
    curve = 2;
    Alpha12 = VLE; 
    xVLE = 0:0.001:1;
    yVLE = zeros(length(xVLE));
    
    for kk = 1:length(xVLE);
        yVLE(kk) = (Alpha12*xVLE(kk))/(1 + xVLE(kk)*(Alpha12 - 1));
    end
end

%%%% Calculate D and B Flow Rates %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B = F * (Xd - Z)/(Xd - Xb);
D = F - B;

%%%% Calculate Operating Lines Based on Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% If Reflux and Feed Condition Specified %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if M == 0;
    
%%%% Calculate Top Operating Line (TOL) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    slopeTOL     = reflux/(reflux + 1);
    interceptTOL = (1/(reflux + 1)) * Xd;

%%%% Calculate Feed Line (FOL) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    slopeFOL     = q/(q - 1);
    interceptFOL = (1/(1 - q)) * Z;

%%%% Find TOL and FOL Intersection %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Xx = (interceptFOL - interceptTOL) /(slopeTOL - slopeFOL);
    
    if q1 == 1
        Xx = Z;
    end
    
    Yy = slopeTOL * Xx + interceptTOL;

%%%% Write TOL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    xTOL = Xx:0.001:Xd;
    yTOL = slopeTOL * xTOL + interceptTOL;

%%%% Write FOL Depending on Feed Conditon %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if Xx < Z
        xFOL = Xx:0.001:Z;
        yFOL = slopeFOL * xFOL + interceptFOL;
    elseif Xx > Z
        xFOL = Z:0.001:Xx;
        yFOL = slopeFOL * xFOL + interceptFOL;
    else
        xFOL = [Xx Xx];
        yFOL = [Z Yy];
    end
 
%%%% Calculate Bottom Operating Line %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    slopeBOL     = (Yy - Xb)/(Xx - Xb);
    interceptBOL = Xb - Xb * slopeBOL;
    xBOL         = Xb:0.001:Xx;
    yBOL         = slopeBOL * xBOL + interceptBOL;
    boilup       = 1/(slopeBOL - 1);

%%%% If reflux and boilup specified %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif M == 1
    
%%%% Calculate Top Operating Line (TOL) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    slopeTOL     = reflux / (reflux + 1);
    interceptTOL = (1 / (reflux + 1)) * Xd;
    
%%%% Calculate Bottom Operating Line (BOL) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    slopeBOL     = (boilup + 1) / boilup;
    interceptBOL = - Xb / boilup;

%%%% Find TOL and BOL Intersection %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Xx = (interceptBOL - interceptTOL) /(slopeTOL - slopeBOL);
    Yy = slopeTOL * Xx + interceptTOL;

%%%% Write TOL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    xTOL = Xx:0.001:Xd;
    yTOL = slopeTOL * xTOL + interceptTOL;

%%%% Write FOL Depending on Feed Conditon %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if Xx < Z
        xFOL         = Xx:0.001:Z;
        slopeFOL     = (Yy - Z) / (Xx - Z);
        interceptFOL = Yy - slopeFOL * Xx;
        q            = slopeFOL / (slopeFOL - 1);
        yFOL         = slopeFOL * xFOL + interceptFOL;
    elseif Xx > Z
        xFOL = Z:0.001:Xx;
        slopeFOL     = (Yy - Z) / (Xx - Z);
        interceptFOL = Yy - slopeFOL * Xx;
        q            = slopeFOL / (slopeFOL - 1);
        yFOL         = slopeFOL * xFOL + interceptFOL;
    else
        xFOL = [Xx Xx];
        yFOL = [Z Yy];  
        q    = 1;
    end
   
%%%% Write BOL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    xBOL = Xb:0.01:Xx;
    yBOL = slopeBOL * xBOL + interceptBOL;

%%%% If boilup and feed condition specified %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif M == 2

%%%% Calculate Feed Line (FOL) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    slopeFOL     = q/(q - 1);
    interceptFOL = (1/(1 - q)) * Z;
    
%%%% Calculate Bottom Operating Line (BOL) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    slopeBOL     = (boilup + 1) / boilup;
    interceptBOL = - Xb / boilup;

%%%% Find BOL and FOL Intersection %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Xx = (interceptFOL - interceptBOL) /(slopeBOL - slopeFOL);
    
    if q1 == 1
        Xx = Z;
    end 
    
    Yy = slopeBOL * Xx + interceptBOL;

%%%% Write BOL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    xBOL = Xb:0.001:Xx;
    yBOL = slopeBOL * xBOL + interceptBOL;

%%%% Write FOL Depending on Feed Conditon %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if Xx < Z
        xFOL = Xx:0.001:Z;
        yFOL = slopeFOL * xFOL + interceptFOL;
    elseif Xx > Z
        xFOL = Z:0.001:Xx;
        yFOL = slopeFOL * xFOL + interceptFOL;
    else
        xFOL = [Xx Xx];
        yFOL = [Z Yy];
    end

%%%% Calculate Top Operating Line %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    slopeTOL     = (Yy - Xd)/(Xx - Xd);
    interceptTOL = Xd - Xd * slopeTOL;
    xTOL         = Xx:0.001:Xd;
    yTOL         = slopeTOL * xTOL + interceptTOL;
    reflux       = slopeTOL /(1 - slopeTOL);

end

%%%% Clear plot and set up axes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold off
hold on
title('McCabe-Thiele Diagram at Specified Parameters')
xlabel('Liquid Phase Mole Fraction [x1]')
ylabel('Vapor Phase Mole Fraction [y1]')
axis([0 1 0 1]);

%%%% Setup x-y Line %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xyx = 0:0.001:1;
xyy = 0:0.001:1;

%%%% Plot Lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(xyx, xyy, 'k')
plot(xVLE, yVLE, 'k')
plot(xTOL, yTOL, 'g')
plot(xBOL, yBOL, 'r')
plot(xFOL, yFOL, 'c')


% Create Stage Steps
if stages == 1
    pause
    % Initialize variables
    X = Xd;
    Y = Xd;
    zz = 1;
    pinch = 0;
    
    % Determine "step-over" location
    step_over = Xx;

    % Step to TOL
    while X(zz) >= step_over && pinch < 40
        zz    = zz + 1;
        pinch = pinch + 1;
        
        % If Constant Relative Volatility        
        if curve ~= 1;       
            X(zz) = Y(zz - 1) / (Alpha12*(1 - Y(zz-1)) + Y(zz-1));
            
        % If VLE Curve Fit    
        else                 
            EQ = VLE_EQ - Y(zz - 1);
            S = solve(EQ, K);
            for jj = 1:length(S)
                if S(jj) > 0 && S(jj) < 1 && isreal(S(jj))
                    X(zz) = double(S(jj));
                end
            end
        end
        
        Y(zz) = slopeTOL * X(zz - 1) + interceptTOL;        
    end
    
    % Display error if pinch point exists
    if pinch >= 40
        error('Pinch point or impossible parameters');
    end
    
    % Determine Feed Location
    FeedLoc = length(X)/2;
    
    % Step to BOL
    while Xb <= X(zz)
        zz    = zz + 1;
        
        % If Constant Relative Volatility
        if curve ~= 1;       
            X(zz) = Y(zz - 1) / (Alpha12*(1 - Y(zz-1)) + Y(zz-1));
        
        % If VLE Curve Fit    
        else                 
            EQ = VLE_EQ - Y(zz - 1);
            S = solve(EQ, K);
            for jj = 1:length(S)
                if S(jj) > 0 && S(jj) < 1 && isreal(S(jj))
                    X(zz) = double(S(jj));
                end
            end
        end
        
        Y(zz) = slopeBOL * X(zz - 1) + interceptBOL;
    end
    
    % Make Final Step to y-x Line
    X(length(X)+1) = X(length(X));
    Y(length(Y)+1) = X(length(X));
        
    % Calculate Partial Stage Requirement
    Xlength = length(X);
    parStage = (X(Xlength-2) - Xb)/(X(Xlength - 2) - X(Xlength - 1));
    
    % Calculate Total Number of Stages
    stage_tot  = floor(length(X)/2);
    if parStage >= 0.01
        stage_theo = stage_tot - 1 + parStage;
    else 
        stage_theo = stage_tot;
    end
    
    % Calculate Distillate and Bottoms Concentrations
    Xbottoms    = X(length(X));
    Xdistillate = Xd;    
    
    % Plot Steps
    for vv = 2:(length(X))
        stepx = [X(vv-1) X(vv)];
        stepy = [Y(vv-1) Y(vv)];
        plot(stepx, stepy, 'b')
    end
    
    % Calculate Minimum Number of EQ Trays at Infinite Reflux
    Xmin = Xd;
    Ymin = Xd;
    ww = 1;
    while Xmin(ww) >= 1.01*Xb 
        ww = ww + 1;
              
        % If Constant Relative Volatility 
        if curve ~= 1;       
            Xmin(ww) = Ymin(ww - 1) / (Alpha12*(1 - Ymin(ww-1)) + Ymin(ww-1));
            
        % If VLE Curve Fit    
        else                 
            EQ = VLE_EQ - Ymin(ww - 1);
            S = solve(EQ, K);
            for jj = 1:length(S)
                if S(jj) > 0 && S(jj) < 1 && isreal(S(jj))
                    Xmin(ww) = double(S(jj));
                end
            end
        end
        
        Ymin(ww) = Xmin(ww - 1);       
    end
    Nmin = length(Xmin)/2;
    
    % Calculate Minimum Theoretical Stages
    par_stage_min = (Xmin(length(Xmin)-1)-Xb) / (Xmin(length(Xmin)-1)-Xmin(length(Xmin)));
    if par_stage_min >= 0.01
        Nmin_theo = Nmin - 1 + par_stage_min;
    else 
        Nmin_theo = Nmin;
    end
       
end

% Find Minumum Reflux Ratio for Specified Feed Condtion "Pinch"
% Find pinch point if q = 1
if q1 == 1
    pinch_x = Z;
    
    if curve == 2
        pinch_y = (Alpha12*pinch_x)/(1+pinch_x*(Alpha12 - 1));
    else
        pinch_y = p(1)*pinch_x^10 + p(2)*pinch_x^9 + p(3)*pinch_x^8 + p(4)*pinch_x^7 + p(5)*pinch_x^6 + p(6)*pinch_x^5 + p(7)*pinch_x^4 + p(8)*pinch_x^3 + p(9)*pinch_x^2 + p(10)*pinch_x+ p(11);
    end
% Find Pinch Point if q ~= 1   
% If Constant Relative Volatility 

elseif q1 ~= 1 && curve == 2
    syms K2
    slopeFOLmin = (Yy - Z)/(Xx - Z);
    interceptFOLmin = Z - slopeFOLmin * Z;
    PinchEQ = (Alpha12*K2)/(1+K2*(Alpha12 - 1)) - (slopeFOLmin * K2) - interceptFOLmin;
    S1 = solve(PinchEQ, K2);
    
    for jjj = 1:length(S1)
        if S1(jjj) > 0 && S1(jjj) < 1 && isreal(S1(jjj))
           pinch_x = double(S1(jjj));
        end
    end
    
    pinch_y = slopeFOLmin * pinch_x + interceptFOLmin;
                
% If VLE Curve Fit

else
    syms K2
    slopeFOLmin = (Yy - Z)/(Xx - Z);
    interceptFOLmin = Z - slopeFOLmin * Z;
    VLE_EQ2 = p(1)*K2.^10 + p(2)*K2.^9 + p(3)*K2.^8 + p(4)*K2.^7 + p(5)*K2.^6 + p(6)*K2.^5 + p(7)*K2.^4 + p(8)*K2.^3 + p(9)*K2.^2 + p(10)*K2+ p(11)-slopeFOLmin*K2 - interceptFOLmin;
    S2 = solve(VLE_EQ2, K2);
    
    for jjj = 1:length(S2)
        if S2(jjj) > 0 && S2(jjj) < 1 && isreal(S2(jjj))
           pinch_x = double(S2(jjj));
        end
    end
    
    pinch_y = slopeFOLmin * pinch_x + interceptFOLmin;
end

% Prevent Negative Reflux and Negative Boilup
if pinch_x <= Xb
    pinch_x = Xb;
end

if pinch_y >= Xd
    pinch_y = Xd;
end

    


% Calculate Minimum Reflux
slopeTOLmin = (pinch_y - Xd) / (pinch_x - Xd);
min_reflux = slopeTOLmin / (1 - slopeTOLmin);

% Calculate Minimum Boilup
slopeBOLmin = (pinch_y - Xb) / (pinch_x - Xb);
min_boilup = 1 / (slopeBOLmin - 1);


% Print Column Specifications %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n')
fprintf('==============================================================\n')
fprintf('==============================================================\n')

% Print Specs from Stage-by-Stage if used
if stages == 1
    fprintf('Actual Distillate Concentration  : %6.2f\n', Xdistillate)
    fprintf('Actual Bottoms Concentration     : %6.2f\n', Xbottoms)
    fprintf('==============================================================\n')
    fprintf('Number of Equilibrium Stages     : %6.0f\n', stage_tot)
    fprintf('Number of Theoretical Stages     : %6.2f\n', stage_theo)
    fprintf('==============================================================\n')
    fprintf('Minimum Number of Stages [Nmin]  : %6.0f\n', Nmin)
    fprintf('Theoretical Minimum Stages       : %6.2f\n', Nmin_theo)
    fprintf('==============================================================\n')
    fprintf('Minimum Reflux Ratio [L/D]min    : %6.2f\n', min_reflux)
    fprintf('Minimum Boilup Ratio [Vbar/B]min : %6.2f\n', min_boilup)
    fprintf('==============================================================\n')    
    fprintf('Feed Stage                       : %6.0f\n', FeedLoc)
    fprintf('==============================================================\n')
end

% Continue to print column specs
fprintf('Feed Condition [q]               : %6.2f\n', q)
fprintf('==============================================================\n')
fprintf('External Reflux Ratio [L/D]      : %6.2f\n', reflux)
fprintf('External Boilup Ratio [Vbar/B]   : %6.2f\n', boilup)
fprintf('==============================================================\n')
fprintf('Overhead Flow Rate [D]           : %6.2f kmol/hr \n', D)
fprintf('Bottoms  Flow Rate [B]           : %6.2f kmol/hr \n', B)
fprintf('==============================================================\n')
fprintf('==============================================================\n')
fprintf('\n\n\n\n\n\n\n')
if stages ~= 1 
    fprintf('\n\n\n\n\n\n\n\n\n\n\n\n\n')
end

% Plot Nmin Steps
if minstep == 1
    pause
    g = figure;
    set(g,'name','Infinte Reflux','numbertitle','off')
    hold on
    title('McCabe-Thiele Diagram at Infinte Reflux')
    xlabel('Liquid Phase Mole Fraction [x1]')
    ylabel('Vapor  Phase Mole Fraction [y1]')
    axis([0 1 0 1]);
    plot(xyx, xyy, 'k')
    plot(xVLE, yVLE, 'k')
    Xmin(length(Xmin)+1) = Xmin(length(Xmin));
    Ymin(length(Ymin)+1) = Xmin(length(Xmin));
    plot(Xd, Xd, 'mo')
    plot(Xb, Xb, 'mo')
    
    for vv = 2:(length(Xmin))
        stepx_min = [Xmin(vv-1) Xmin(vv)];
        stepy_min = [Ymin(vv-1) Ymin(vv)];
        plot(stepx_min, stepy_min, 'b')
    end
    
end

% Plot Minimum Reflux Diagram
if minstep == 1 && maxstep == 1
    pause
    q = figure;
    set(q,'name','Minimum Reflux','numbertitle','off')
    hold on
    title('McCabe-Thiele Diagram at Minimum Reflux for Given Feed')
    xlabel('Liquid Phase Mole Fraction [x1]')
    ylabel('Vapor  Phase Mole Fraction [y1]')
    axis([0 1 0 1]);
    plot(xyx, xyy, 'k')
    plot(xVLE, yVLE, 'k')
    plot(Xd, Xd, 'mo')
    plot(Xb, Xb, 'mo')
    
    % Create TOL
    minTOLx = [pinch_x Xd];
    minTOLy = [pinch_y Xd];
    
    % Create BOL
    minBOLx = [Xb pinch_x];
    minBOLy = [Xb pinch_y];
    
    % Create FOL
    if pinch_x < Z
        minFOLx = [pinch_x Z];
        minFOLy = [pinch_y Z];
        
    elseif pinch_x > Z
        minFOLx = [Z pinch_x];
        minFOLy = [Z pinch_y];
        
    else
        minFOLx = [Z Z];
        minFOLy = [Z pinch_y];
    end
    
    % Plot Operating Lines
    plot(minTOLx , minTOLy, 'g')
    plot(minBOLx , minBOLy, 'r')
    plot(minFOLx , minFOLy, 'c')
    
end

pause
% Make Column Diagram
if column == 1;
    c = figure;
    set(c, 'name', 'Column Diagram', 'numbertitle', 'off')
    hold on
    title('Column at Specs [assuming total condenser and partial reboiler')
    axis([0 100 0 100])
    
    col_base = 20;
    col_height = col_base + 60;
    col_width = col_base + 20;
    trays = stage_tot;
    feed_tray = FeedLoc;
    tray_height = (col_height - col_base) / (trays - 1);
    tray_point = col_height - 0.5*tray_height;
    bottoms_point = col_base + 0.5*tray_height;
    
    % Column Frame
    % Left Wall
    LWx = [col_base col_base];
    LWy = [col_base col_height];
    
    % Feed Location
    
    % Right Wall
    RWx = [col_width col_width];
    RWy = [col_base col_height];
    
    % Top Wall
    TWx = [col_base col_width];
    TWy = [col_height col_height];
    
    % Bottom Wall
    BWx = [col_base col_width];
    BWy = [col_base col_base];
    
    % Total Condenser
    up_line1 = 0.8*(col_width-col_base) + col_base;
    up_line2 = col_height + 10;
    upline_x = [up_line1 up_line1];
    upline_y = [col_height up_line2];
    rightline1 = [up_line1 (up_line1 + 15)];
    rightline2 = [up_line2 up_line2];
    r = 2.5;
    condenser_centerx = rightline1(2) + r;
    condenser_centery = up_line2;
    slash_x = [rightline1(2) (rightline1(2)+2*r)];
    slash_y = [(rightline2(2)-1.5*r) (rightline2(2)+1.5*r)];
    arrow_left_x = [slash_x(1) slash_x(1)];
    arrow_right_x = [slash_x(1) (slash_x(1)+2)];
    arrow_left_y = [slash_y(1) (slash_y(1)+2)];
    arrow_right_y = [slash_y(1) slash_y(1)];
    angles = 0:0.001:2*pi;
    xp=r*cos(angles);
    yp=r*sin(angles);
    downlinex = [condenser_centerx condenser_centerx];
    downliney = [(condenser_centery - r) tray_point];
    horizontal_linex = [col_width (col_width + 20)];
    horizontal_liney = [tray_point tray_point];
    arrow3x = [col_width col_width+3];
    arrow4x = [(col_width + 20) (col_width +17)];
    arrow3y = [tray_point (tray_point+2)];
    arrow3y2 = [tray_point (tray_point-2)];
    
    plot(arrow_left_x, arrow_left_y,'k')
    plot(arrow_right_x,arrow_right_y,'k')
    plot(arrow3x, arrow3y,'k')
    plot(arrow4x, arrow3y, 'k')
    plot(arrow4x, arrow3y2, 'k')
    plot(arrow3x, arrow3y2,'k')
    plot(horizontal_linex, horizontal_liney,'k')
    plot(downlinex, downliney, 'k')
    plot(condenser_centerx+xp,condenser_centery+yp, 'k');
    plot(slash_x, slash_y, 'k');
    plot(rightline1, rightline2, 'k')
    plot(upline_x, upline_y, 'k')
    
    % Partial Reboiler
    up_line3 = 0.8*(col_width-col_base) + col_base;
    up_line4 = col_base - 10;
    upline_x1 = [up_line3 up_line3];
    upline_y1 = [col_base up_line4];
    rightline3 = [up_line3 (up_line3 + 15)];
    rightline4 = [up_line4 up_line4];
    reboiler_centerx = rightline3(2) + r;
    reboiler_centery = up_line4;
    upslash_x = [rightline3(2) (rightline3(2)+2*r)];
    upslash_y = [(rightline4(2)-1.5*r) (rightline4(2)+1.5*r)];
    arrow_left_x = [upslash_x(1) upslash_x(1)];
    arrow_right_x = [upslash_x(1) (upslash_x(1)+2)];
    arrow_left_y = [upslash_y(1) (upslash_y(1)+2)];
    arrow_right_y = [upslash_y(1) upslash_y(1)];
    xp1=r*cos(angles);
    yp1=r*sin(angles);
    downlinex = [reboiler_centerx reboiler_centerx];
    downliney = [(reboiler_centery + r) bottoms_point];
    horizontal_linex = [col_width reboiler_centerx];
    horizontal_liney = [bottoms_point bottoms_point];
    arrow3x  = [col_width col_width+3];
    arrow5x  = [(reboiler_centerx + r) (reboiler_centerx + 10)];
    arrow5y  = [reboiler_centery reboiler_centery];
    arrow3y  = [bottoms_point (bottoms_point+2)];
    arrow3y2 = [bottoms_point (bottoms_point-2)];
    arrow7x  = [reboiler_centerx+r-1.5 reboiler_centerx+r];
    arrow7x2 = [arrow7x(2) arrow7x(2)];
    arrow7y1 = [upslash_y(2) upslash_y(2)];
    arrow7y2 = [upslash_y(2) (upslash_y(2)-2.5)];
    arrow6x  = [arrow5x(2)-2 arrow5x(2)];
    arrow6y1 = [arrow5y(1)+2 arrow5y(1)];
    arrow6y2 = [arrow5y(1)-2 arrow5y(1)];
    arrow8x  = [(reboiler_centerx-r-2) (reboiler_centerx-r)];
    arrow8y1 = [(reboiler_centery+2) reboiler_centery];
    arrow8y2 = [(reboiler_centery-2) reboiler_centery];
    arrow9x  = [(condenser_centerx-r-2) (condenser_centerx-r)];
    arrow9y1 = [(condenser_centery+2) condenser_centery];
    arrow9y2 = [(condenser_centery-2) condenser_centery];
    
    plot(arrow8x, arrow8y1,'k')
    plot(arrow8x, arrow8y2,'k')
    plot(arrow9x, arrow9y1,'k')
    plot(arrow9x, arrow9y2,'k')
    plot(arrow6x, arrow6y1,'k')
    plot(arrow6x, arrow6y2,'k')
    plot(arrow7x, arrow7y1,'k')
    plot(arrow7x2, arrow7y2,'k')
    plot(arrow5x, arrow5y,'k')
    plot(arrow3x, arrow3y2,'k')
    plot(arrow3x, arrow3y,'k')
    plot(horizontal_linex, horizontal_liney,'k')
    plot(downlinex, downliney, 'k')
    plot(reboiler_centerx+xp1,reboiler_centery+yp1, 'k');
    plot(upslash_x, upslash_y, 'k');
    plot(rightline3, rightline4, 'k')
    plot(upline_x1, upline_y1, 'k')
    
    % Plot Stage Trays
    tray = col_height;
    tray_x = [col_base col_width];
    note_position = 0.5*(col_width - col_base);
    
    
    for nn = 1:(trays-1)
        tray = tray - tray_height;
        tray_y = [tray tray];
        plot(tray_x, tray_y, 'b') ;
        if nn == feed_tray && nn ~= trays
            feed_x = [10 20];
            feed_y = [(0.5*tray_height + tray) (0.5*tray_height + tray)];
            plot(feed_x, feed_y, 'r')
            arrow_x = [17 20];
            arrow_y_top = [(2+feed_y(1)) feed_y(1)];
            arrow_y_bot = [(feed_y(1)-2) feed_y(1)];
            plot(arrow_x, arrow_y_top, 'r')
            plot(arrow_x, arrow_y_bot, 'r')
        end
        
    end
    if feed_tray == trays
        feed_x = [10 reboiler_centerx];
        feed_y = [3 3];
        feedupx = [reboiler_centerx reboiler_centerx];
        feedupy = [3 reboiler_centery-r];
        feedup_x1 = [feedupx(1)-1.5 feedupx(1)];
        feedup_x2 = [feedupx(1) feedupx(1)+1.5];
        feedup_y1 = [feedupy(2)-1.5 feedupy(2)];
        feedup_y2 = [feedupy(2) feedupy(2)-1.5];
        downlinex = [10 10];
        downliney = [30 3];
        rlinex = [5 10];
        rliney = [30 30];
        
        plot(rlinex, rliney, 'r')
        plot(downlinex, downliney, 'r')
        plot(feed_x, feed_y,'r')
        plot(feedupx, feedupy, 'r')
        plot(feedup_x1,feedup_y1,'r')
        plot(feedup_x2,feedup_y2,'r')
    end
    
    
    %Plot Frame
    plot(LWx, LWy, 'k')
    plot(RWx, RWy, 'k')
    plot(TWx, TWy, 'k')
    plot(BWx, BWy, 'k')
    
    
end

end

% VLE     = [0 0; 0.01 0.03; 0.02 0.134; 0.04 0.23; 0.06 0.304; 0.08 0.365; 0.1 0.418; 0.15 0.517; 0.2 0.579; 0.3 0.665; 0.4 0.729; 0.5 0.779; 0.6 0.825; 0.7 0.87; 0.8 0.915; 0.9 0.958; 0.95 0.979; 1 1];




Contact us