function [sys,x0] = colas(t,x,u,flag)
%
%  Simulink interface to colamod.m     
%
% Inputs:    t    - time in [min].
%            X    - State, the first 41 states are compositions of light
%                   component A with reboiler/bottom stage as X(1) and 
%                   condenser as X(41). State X(42)is holdup in reboiler/
%                   bottom stage and X(82) is hold-up in condenser. 
%            U(1) - reflux L,
%            U(2) - boilup V,
%            U(3) - top or distillate product flow D,
%            U(4) - bottom product flow B,
%            U(5) - feed rate F,
%            U(6) - feed composition, zF.
%            U(7) - liquid feed fraction, qF.
%
% Outputs:   sys and x0 as described in the SIMULINK manual.
%            when flag is 0 sys contains sizes and x0 contains 
%            initial condition. 
%            when flag is 1, sys contains the derivatives,
%            and when flag is 3 sys contains outputs; 
%            y(1)    - top composition,
%            y(2)    - bottom composition,
%            y(3)    - condenser holdup,
%            y(4)    - reboiler holdup,
%            y(5:45) - tray composition, y(5) is reboiler y(45) is top.

NT = 41;

Tb_ben=353.3;
Tb_tol=383.8;
Tb_xyl= 411.5;




if abs(flag) == 1
  % Return state derivatives.
  sys = colamod(t,x,u);
elseif abs(flag) == 3
  % Return system outputs.
  sys(1,1) = x(41);         % Top composition.
  sys(2,1) = x(205);          % Bottom composition.
  sys(3,1) = 1 - x(124)-x(165);       % Holdup in condenser.
  sys(4,1) = x(123);       % Holdup in condenser.
  sys(5,1) = x(83);       % Holdup in condenser.
  sys(6,1) = x(246);       % Holdup in condenser.
  sys(7,1) = x(206);       % Holdup in reboiler.
  
    sys(8:41+7,1)= Tb_ben*x(1:41)+(x(42:82)*Tb_tol) + ((1-x(1:41)-x(42:82))*Tb_xyl) ; %Temperatures
  sys(49:41+48,1) = x(1:41) ;
  sys(90:41+89,1) = x(42:82) ;
elseif flag == 0
  % Ininitalize the system
LT1 = 3.28852 ;
LT2 = 1.66579 ;
VB1 = 3.86453 ;
VB2 = 1.89747 ;
B1 = 0.81398  ;
B2 = 0.82024 ;
D1 = 0.57602 ;
D2 = 0.57976 ;
qF = 1;
zF = [0.4 0.2 0.4];
F = 1;

U(1)=LT1 ;
U(2)=VB1 ;
U(3)=LT2 ; 
U(4)=VB2 ; 
U(5)=D1 ; 
U(6)=B1 ;
U(7)=D2 ; 
U(8)=B2 ;
U(9)= F ;
U(10)= qF;



x0 =[0.01079 0.01720 0.02542 0.03556 0.04770 0.06191 0.07828 0.09690 0.11780 0.14096 0.16623 0.19337 0.22197 0.25152 0.28142 0.31102 0.33971 0.36693 0.39224 0.41534 0.43607 0.47736 0.51535 0.54977 0.58104 0.60990 0.63710 0.66329 0.68891 0.71422 0.73930 0.76413 0.78859 0.81250 0.83568 0.85791 0.87903 0.89888 0.91736 0.93441 0.95000 ;
0.30616 0.37920 0.44259 0.49313 0.52982 0.55339 0.56549 0.56807 0.56296 0.55175 0.53578 0.51616 0.49388 0.46986 0.44492 0.41983 0.39528 0.37184 0.34996 0.32994 0.31196 0.33665 0.35184 0.35795 0.35624 0.34821 0.33530 0.31874 0.29951 0.27839 0.25602 0.23293 0.20957 0.18636 0.16363 0.14167 0.12072 0.10097 0.08255 0.06554 0.04998 ]';
x1(1:41)=x0(1:41,1)';
x1(42:82)=x0(1:41,2)';
x01 =[0.00001 0.00001 0.00002 0.00003 0.00004 0.00006 0.00009 0.00013 0.00018 0.00025 0.00035 0.00047 0.00064 0.00085 0.00114 0.00151 0.00199 0.00262 0.00343 0.00448 0.00583 0.00600 0.00620 0.00644 0.00672 0.00704 0.00740 0.00778 0.00819 0.00862 0.00908 0.00958 0.01017 0.01089 0.01186 0.01320 0.01514 0.01799 0.02222 0.02852 0.03790 ;
0.04999 0.06772 0.08692 0.10735 0.12867 0.15049 0.17237 0.19387 0.21457 0.23412 0.25224 0.26875 0.28355 0.29662 0.30800 0.31778 0.32606 0.33297 0.33863 0.34316 0.34664 0.37017 0.39888 0.43316 0.47302 0.51797 0.56692 0.61826 0.67003 0.72018 0.76691 0.80886 0.84526 0.87585 0.90078 0.92043 0.93523 0.94557 0.95162 0.95327 0.95000 ]';
x2(1:41)=x01(1:41,1)';
x2(42:82)=x01(1:41,2)';


x1(83:123) =[0.53140 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.54200 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.53468 0.50760 ]';
x2(83:123) =[0.50823 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.44124 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.43090 0.47317 ]';

x1 = x1';
x2 = x2';

x0 = [x1;x2] ;
 
  sys = [6*NT, 0, 130, 10, 0, 0];
else
  sys = [];  
end

