%**************************************************************
%            A PROGRAM FOR RUNNING THE HP COLUMN
%
% - calls on colamod1.m the model for the HP column
% - The no. of stages, feed data must be changed for a 
%   new case
%
% H.K.Engelien, 23.10.2001
%**************************************************************

function xprime = col_pressure(t,X,options,u0,s0)

global Ptop_HP Y1 time Flow 

NT = 16;                         % No. stages in HP column
NF = 4;                          % Feed stage (counted from btm)

X1 = X(1:3*NT);                  % States in HP column

% Feed
F_HP  = 4320;              % Feedrate (kmol/h)
zF_HP = 0.73;              % Feed composition light component
qF_HP = 1;                 % Feed liq. fraction
TF    = 394.65;             % Feed temperature (K)

% Get actual values
MB_HP = X1(NT+1);          % Actual reboiler holdup (kmol)
MD_HP = X1(2*NT);          % Actual condenser holdup (kmol)
xD_HP = X1(NT);            % Concentration in distillate
xB_HP = X1(1);             % Concentration in btm product
TD_HP = X1(3*NT);          % Temperature in condenser (K)
TB_HP = X1(2*NT+1);
T(NT) = X1(3*NT);
x(NT) = X1(NT);

% Calculating top pressure for cooling control:

% Wilson equation coefficients :
% (from Gmehling & Onken VLE Data Collection, p.51)
A11 = 1;
A12 = -102.7135;           % Reference pressure : 5.066 bar
A21 = 695.4902;            % Reference pressure : 5.066 bar     
A22 = 1;
V_water    = 18.07;        % Molar vol. of pure liq. component
V_methanol = 40.73;        % Molar vol. of pure liq. component 

% Antoine Coefficients
A1 = 18.5875; B1 = 3626.55; C1 = -34.29;
A2 = 18.3036; B2 = 3816.44; C2 = -46.13;

% Activity coefficients from Wilson equation
i=NT;
R = 1.98721;                  % Universal gas constant (cal/kmol*K)
lambda_12(i) = V_water/V_methanol*exp(-(A12-A11)./(R.*T(i)));
lambda_21(i) = V_methanol/V_water*exp(-(A21-A22)./(R.*T(i)));
gammaL(i)    = exp(-log(x(i)+lambda_12(i).*(1-x(i)))...
               +(1-x(i)).*(lambda_12(i)./(x(i)+lambda_12(i).*(1-x(i)))...
               -lambda_21(i)./(lambda_21(i).*x(i)+(1-x(i)))));
gammaH(i)    = exp(-log((1-x(i))+lambda_21(i).*x(i))...
               -x(i).*(lambda_12(i)./(x(i)+lambda_12(i).*(1-x(i)))...
               -lambda_21(i)./(lambda_21(i).*x(i)+(1-x(i)))));

% Vapour pressure fra Antoine equation
i=NT;
pL(i) = exp(A1 - B1./(T(i)+C1))*133.32;            % Vap. pressure light (Pa) 
pH(i) = exp(A2 - B2./(T(i)+C2))*133.32;            % Vap. pressure heavy (Pa)
P_top =  x(i).*pL(i).*gammaL(i) + (1-x(i)).*pH(i).*gammaH(i);          % Stage pres.(Raoult's law)

% PI-controllers for Control of HP Column:
%-------------------------------------------------------------
 
D0_HP  = u0(1);                    % Nominal distillate flow (kmol/h) 
B0_HP  = u0(2);                    % Nominal bottom flow (kmol/h)
VD0_HP = u0(3);                    % Nominal flow condensated (kmol/h)
LT0_HP = u0(4);                    % Nominal reflux (kmol/h)
VB0_HP = u0(5);                    % Nominal boilup (kmol/h)

% Top (condenser) holdup control
KcD  = -60;                       
MDs  = s0(1);
D_HP = D0_HP + (MDs-MD_HP)*KcD;  

% Bottom (reboiler) holdup control
KcB  = -60;                                  
MBs  = s0(2);         
B_HP = B0_HP + (MBs-MB_HP)*KcB;    % Uses btm. flow to control holdup

% Distillate composition control
KcXD    = 500000;
TauI_xD = 0.067;
xDs_HP  = s0(4);
xD_HP_error  = (xDs_HP-xD_HP); Interr_xD = X(3*NT+1);
LT_HP   = LT0_HP + KcXD*xD_HP_error+KcXD/TauI_xD*Interr_xD;  

% Btm composition control
KcXB    = -50000;
TauI_xB = 0.055;      
xBs_HP  = s0(3);
xB_HP_error = (xBs_HP-xB_HP);Interr_xB = X(3*NT+2);
VB_HP   = VB0_HP + KcXB*xB_HP_error + KcXB/TauI_xB*Interr_xB; 

% Pressure control (top)
KcP     = -0.028;
TauI_P  = 0.067;
Ps_HP   = s0(5);           
P_HP_error = (Ps_HP-P_top); Interr_P = X(3*NT+3);
VD_HP   = VD0_HP + KcP*P_HP_error + KcP/TauI_P*Interr_P;

%---------------------------------------------------------------------

% Store all inputs and disturbances (goes into colamod1.m)
U1(1)=LT_HP;U1(2)=VB_HP;U1(3)=D_HP;U1(4)=B_HP;U1(5)=F_HP;
U1(6)=zF_HP;U1(7)=qF_HP;U1(8)=VD_HP;U1(9)=TF;

xprime1 = colamod1(t,X1,U1);
t
time = [time t];                     % Time vector

Ptop_HP = [Ptop_HP; P_top];           % Top pressure response 

Y1 =[D_HP B_HP VD_HP LT_HP VB_HP xD_HP xB_HP];

Flow = [Flow; D_HP B_HP VD_HP LT_HP VB_HP];


IntegralError= [xD_HP_error;xB_HP_error;P_HP_error];  %  Integral error
xprime = [xprime1; IntegralError];


