%**************************************************************
% ******************************************************************
%   A PROGRAM FOR RUNNING A COLUMN WITH PRESSURE VARIATIONS
%
% - dynamic simulation of a distillation column with varying pressure
% - calls on col_pressure.m 
% - inputs:  X - a vector with the states of the column
%            (composition, holdup and temperature)
%            u0 - a vector with the nominal values for the 
%            contol loops in the column
%            s0 - a vector with the setpoints
%
% H.K.Engelien, 23.10.2001
%*******************************************************************

global Ptop_HP Y1 time P_HP Flow  Hvap

% Number of stages in column
NT = 16;        

Ptop_HP = [ ];time = [ ];Flow = [];

% Get input data: nominal operating values, set points
load steady_state
u0 = [D_HP B_HP VD_HP LT_HP VB_HP];    % nominal values for control loops
s0 = [X(32) X(17) 0.01 0.99 P_top_HP]; % setpoints for control loops
Xss= [X;0;0;0];

options = [];
T_end = 10;
[t,x] = ode15s('col_pressure', [0 T_end] ,Xss,options,u0,s0);

t_final = size(x);

X1final = [x(t_final(1),1:3*NT)]';     % Final states HP col

% Display results on screen
X = [X1final];
D_HP= Y1(1), B_HP= Y1(2), VD_HP= Y1(3), LT_HP= Y1(4)
VB_HP= Y1(5), xD_HP= Y1(6),xB_HP= Y1(7)
c=size(Ptop_HP);
P_top_HP = Ptop_HP(c(1))
T_top_HP = x(t_final(1),3*NT)
Q_top_HP = Hvap(NT)*VD_HP*1000/(3600*1000)    % kW
Q_btm_HP = Hvap(1)*VB_HP*1000/(3600*1000)     % kW

%save steady_state X D_HP B_HP VD_HP LT_HP VB_HP xD_HP xB_HP P_top_HP

% Plotting HP column profiles:
%------------------------------------

figure(1)
subplot(2,2,1); plot(x(t_final(1),1:NT),'r-')
title('Composition Profile in HP Column')
subplot(2,2,2); plot(x(t_final(1),NT+1:2*NT),'r-')      
title('Holdup Profile in HP Column')
subplot(2,2,3); plot(x(t_final(1),2*NT+1:3*NT),'r-')      
title('Temperature Profile in HP Column')
subplot(2,2,4); plot(P_HP,'r-')
title('Pressure Profile in HP Columns')
hold;

% Plotting HP column responses over the whole simulation time:
%-------------------------------------------------------------
% Pressure and condensate flow response
figure(2)
subplot(2,1,1); plot(time,Ptop_HP/1e5);hold
plot([0 T_end],[s0(5)/1e5 s0(5)/1e5],'r--');
title('Pressure Response in Top of Column')
ylabel('Pressure (bar)')
subplot(2,1,2); plot(time,Flow(:,3))
title('Condensation (VD) in HP Column')
hold;

% Bottom & top composition response & internal flow response
figure(3)
subplot(2,2,1); plot(t,x(:,1)); hold; plot([0 T_end],[s0(3) s0(3)],'r--');
title('Bottom Composition (LK) Response in HP Column')
subplot(2,2,2); plot(t,x(:,NT)); hold; plot([0 T_end],[s0(4) s0(4)],'r--');
title('Top Composition (LK) Response in HP Column')
subplot(2,2,3); plot(time,Flow(:,5))
title('Boilup Response in HP Column')
subplot(2,2,4); plot(time,Flow(:,4))
title('Reflux Response in HP Column')
hold

% Holdup response  and external flow response
figure(4)
subplot(2,2,1); plot(t,x(:,NT+1));hold; plot([0 T_end],[s0(2) s0(2)],'r--');
title('Btm Holdup Response in HP Column')
subplot(2,2,2); plot(t,x(:,2*NT));hold; plot([0 T_end],[s0(1) s0(1)],'r--');
title('Top Holdup Response in HP Column')
subplot(2,2,3); plot(time,Flow(:,2))
title('Bottom Flowrate Response in HP Column')
subplot(2,2,4); plot(time,Flow(:,1))
title('Distillate Response in HP Column')
hold;




