function [sys,x0,str,ts] = dens_sf(t,x,u,flag,x_0)

switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes(x_0);
  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u);

  %%%%%%%%%%%%%%%%
  % Unused flags %
  %%%%%%%%%%%%%%%%
  case{2,4,5,9}
     sys = [];
     
  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);

end
% end sf_densidade

%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================

function [sys,x0,str,ts]=mdlInitializeSizes(x_0)

% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
sizes = simsizes;
sizes.NumContStates  = 2;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 2;
sizes.NumInputs      = 4;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;   % at least one sample time is needed
sys = simsizes(sizes);

% initialize the initial conditions
x0=x_0;

% str is always an empty matrix
str = [];

% initialize the array of sample times
ts  = [0 0];

% end mdlInitializeSizes

%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================

function sys=mdlDerivatives(t,x,u)

%u(1): loop 1 C3 flow kg/h (MV)
%u(2): loop 1 production rate kg/h (DV)
%u(3): loop 2 C3 flow kg/h (MV)
%u(4): loop 2 production rate kg/h (DV)

%x(1): loop 1 density kg/m3
%x(2): loop 2 density kg/m3

%fixed parameters
roC3 = 400; %kg/m3
roPP = 900; %kg/m3
V = 48;     %m3

sys(1) = x(1)/V/roC3*( (1-roC3/roPP)*u(2) + u(1)*(roC3/x(1)-1) );
sys(2) = x(2)/V/roC3*( (1-roC3/roPP)*u(4) + (u(1)-V*sys(1))*(1-roC3/x(1)) + ...
           (1-roC3/x(2))*(V*sys(1)-u(3)-u(1)) );
         
% end mdlDerivatives

%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================

function sys=mdlOutputs(t,x,u)

sys = x;

% end mdlOutputs
