Code covered by the BSD License

# McCabeThiele

### Vincent Cericola (view profile)

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];