% INIT.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Initializes the Kaibel column

clear all; close all; clc;

global coldata dyndata hydraulic termodata

% COLUMN DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Number of stages
coldata.Nt = [12 12 8 8 8 8 8];
% Height of each section
coldata.sech = [1.33 1.33 0.89 0.89 0.89 0.89 0.89];
% Column diameter
coldata.diameter = 0.1;
% Initial holdup in reboiler and condenser
coldata.M0 = [0.5 0.5 1 0.5 0.5 1 0.5];

% THERMODYNAMIC DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Universal gas constant [cal/molK]
termodata.R = 1.98721;
% Antoine coefficients
termodata.A = [
    8.08097  1582.271  239.726
    8.11220  1592.864  226.184
    8.37895  1788.020  227.438
    7.36366  1305.198  173.427
    ];
% Molar volumes of pure components
termodata.V = [40.73  58.68  75.14  91.97];
% Wilson parameters
termodata.W = [
    0           507.3868    71.4631     -121.3194
    -452.5062   0           102.6706    -8.3924
    154.6476    -40.8745    0           -316.7718
    9543.2503   122.8635    547.4636    0
    ];
% Pressure (atm)
termodata.P = 1.0;
% Heat transfer coefficient (W/m^2 K)
termodata.U = 0;
% Molar enthalpy of vapourisation (J/mol)
termodata.deltaH = 45000;
% Temperature outside column (K)
termodata.T0 = 293;
% Vapour bypass factor
termodata.alpha = 0;

% HYDRAULIC DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hydraulic.k1 = 1/0.063;

% DYNAMIC DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('mat-files/x_nom');
load('mat-files/u_nom');

dyndata.yselect = 1;

dyndata.x_nom = x_nom;
dyndata.u_nom = u_nom;
dyndata.y_nom = genkaibel([],x_nom,u_nom,3,dyndata.yselect);

clear x_nom u_nom

% CONTROLLER DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Controller parameters
% control.X = [
%   Controller gain (Kp)
%   Integrator time (Ti)
%   Maximum controller output (u_max)
%   Minimum controller output (u_min)
%   Timedelay for measurement
%   Windup time constant (Tw)
%   Valve opening time
%   Bias
%   ];
control.onoff = 1;
% Reflux controller
control.L = [
    -0.02005 * control.onoff
    150.2
    0.1508
    -2.8492
    2
    8
    0.1
    0
    ];
% Sidestream 1 controller
control.S1 = [
    0.0058 * control.onoff
    36.8
    0.7506
    -0.2492
    2
    250/10000
    0.1
    0
    ];
% Sidestream 2 controller
control.S2 = [
    0.00475 * control.onoff
    150.7
    0.7503
    -0.2497
    2
    250/10000
    0.1
    0
    ];
% Liquid split controller
control.Rl = [
    -0.0017 * control.onoff
    35.3
    0.7428
    -0.2572
    1
    8
    0.1
    0
    ];
% D controller
control.D = [
    -0.09 * control.onoff
    50
    0.7492
    -0.2508
    2
    25
    0.1
    0];
% B controller
control.B = [
    -0.09 * control.onoff
    50
    0.7497
    -0.2503
    2
    25
    0.1
    0
    ];

if (control.onoff == 1)
    disp('### CONTROLLERS ON ###');
else
    disp('### CONTROLLERS OFF ###');
end