function [sys,x0] = colas_stepvis(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;
Tb1=261.45; %Boling point of iso-butane
Tb2=272.65; %Boling point of n-butane
deltaTb= Tb1-Tb2;

global OP

if abs(flag) == 1
  % Return state derivatives.
  sys = colamod(t,x,u);
elseif abs(flag) == 3
  % Return system outputs.
  sys(1,1)= 1-x(NT);         % Top composition
  %sys(1,1) = -log(x(NT)/(1-x(NT))); % Log transformed top impurity
  sys(2,1)= x(1);          % Bottom composition
  %sys(2,1) = log(x(1)/(1-x(1)));   % Log transformed bottom impurity
  sys(3,1) = x(2*NT);       % Holdup in cond
  sys(4,1) = x(NT+1);       % Holdup in reboiler
  i=5:NT+4; sys(i,1)= deltaTb*x(i-4)+Tb2;    % Temperature at each stage
  sys(NT+5:2*NT+4,1)= x(1:NT);   % Compositions
  sys(2*NT+5:3*NT+4,1)=x(NT+1:NT*2);
  %sys(2*NT+5,1)=deltaTb*x(15)+Tb2;
elseif flag == 0
  % Ininitalize the system
    if OP==1
        load cola_init_1
    elseif OP==2
        load cola_init_2
    elseif OP==3
        load cola_init_3
    elseif OP==4
        load cola_init_4a
    elseif OP==5
        load cola_init_4b
    elseif OP==6
        load cola_init_4c
    elseif OP==7
        load cola_init_5
    elseif OP==8
        load cola_init_6
    elseif OP==9
        load cola_init_7
    else
        load cola_init_8
    end
  x0 = Xinit;
  sys = [2*NT, 0, 3*NT+4, 7, 0, 0];
else
  sys = [];  
end
  
    