function sf_flash(block)
% In order to know more on Level-2 M-file S-function, refer to the MatLab documentation.

setup(block);

function setup(block)

% Register number of dialog parameters
block.NumDialogPrms = 0;

% Register number of ports
nu = 4; ny = 9;
block.NumInputPorts  = nu;
block.NumOutputPorts = ny;

% Setup port properties to be inherited or dynamic
block.SetPreCompInpPortInfoToDynamic;
block.SetPreCompOutPortInfoToDynamic;

for i = 1:nu
    block.InputPort(i).Dimensions = 1;
    block.InputPort(i).DirectFeedthrough = true;
end

for i = 1:ny
    block.OutputPort(i).Dimensions = 1;
end

block.SampleTimes = [0 0];

% Setup Dwork
block.NumContStates = 3;

block.RegBlockMethod('InitializeConditions',    @InitializeConditions);
block.RegBlockMethod('Outputs',                 @Outputs);
block.RegBlockMethod('Derivatives',             @Derivatives);
block.RegBlockMethod('SetInputPortSamplingMode',@SetInpPortFrameData);


function SetInpPortFrameData(block, idx, fd)
  
block.InputPort(idx).SamplingMode = fd;
for i = 1:block.NumOutputPorts
    block.OutputPort(i).SamplingMode = fd;
end


function InitializeConditions(block)

load steadystate % y0
block.ContStates.Data = y0(1:3);


function Outputs(block)

parameters
y = block.ContStates.Data;
for i = 1:block.NumInputPorts
    u(i) = block.InputPort(i).Data;
end

% Disturbances
q  = u(1) + qbias;
Tf = u(2) + Tfbias;
z1 = u(3) + z1bias;

% Manipulable variables
L  = u(4) + Lbias;

% State variables
% II. Extract present value of states
n  = y(1);                               % total holdup in tank (both phases)       [mol]  
n1 = y(2);                               % component 1 holdup in tank (both phases) [mol1]
U  = y(3);                               % total internal energy (both phases)      [J]

% Non-linear model
T = fzero(@myfun,T0,[],n,n1,U,q,Tf,z1,L);
model

% Output variables
yout = [G y1 (1-y1) L x1 (1-x1) p T VL] - youtbias;
for i = 1:block.NumOutputPorts
    block.OutputPort(i).Data = yout(i);
end


function Derivatives(block)

parameters
y = block.ContStates.Data;
for i = 1:block.NumInputPorts
    u(i) = block.InputPort(i).Data;
end

% Disturbances
q  = u(1) + qbias;
Tf = u(2) + Tfbias;
z1 = u(3) + z1bias;

% Manipulable variables
L  = u(4) + Lbias;

% State variables
% II. Extract present value of states
n  = y(1);                               % total holdup in tank (both phases)       [mol]  
n1 = y(2);                               % component 1 holdup in tank (both phases) [mol1]
U  = y(3);                               % total internal energy (both phases)      [J]

% Non-linear model
T = fzero(@myfun,T0,[],n,n1,U,q,Tf,z1,L);
model

block.Derivatives.Data = F;


function g = myfun(T,n,n1,U,q,Tf,z1,L)

parameters

model
g = U + p*Vtot - hL*nL - hG*nG;


