function [sys,x0,str,ts] = separator_2phases(t,x,u,flag,R,rho_gasS,MW_gas,rho_oilS,...
    rho_oil,R_sep,L_sep,n,Tgas,FlashFactor,h_oil_meas,ho0,p0,PAHH,PALL,Oil_LAHH,Oil_LALL,Orientation,Ts,Tag)


% Parameters:
% R                   Gas constant [J/mol K]
% rho_gasS            Gas density at standard conditions [kg/Sm3]
% MW_gas              Molar weight of gas [kg/mol]
% rho_oilS            Oil density at standard conditions [kg/Sm3]
% rho_oil             Oil density [kg/m3]
% rho_waterS          Water density at standard conditions [kg/Sm3]
% rho_water           Water density [kg/m3]
% R_sep               Separator radius [m]
% L_sep               Separator length [m]
% n                   Radii_ratio = n:1. Must be a positive number
% Tgas                Gas temperature [K]
% h_oil_meas          Height of oil outlet [m]
% hw0                 Initial water level [m]
% ho0                 Initial oil level [m]
% p0                  Initial gas pressure [bara]
% PAHH                High gas pressure trip level [bara]
% PALL                High gas pressure trip level [bara]
% Oil_LAHH            High level trip, oil [m]
% Oil_LALL            Low level trip, oil [m]
% Ts                  Sampling time
%
% Inputs:
% 1. Oil rate in [Sm3/s]
% 2. Gas rate out [Sm3/s]
% 3. Oil rate out [m3/s]
%
% Outputs:
% 1. Gas pressure [bara]
% 2. Oil level [m]
% 3. Pressure @ oil output [bara]




switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts] = mdlInitializeSizes(R_sep,L_sep,n,ho0,p0,Orientation);
  %%%%%%%%%%
  % Update %
  %%%%%%%%%%
  case 2,                                                
    sys = mdlUpdate(t,x,u,R,rho_gasS,MW_gas,rho_oilS,rho_oil,R_sep,L_sep,n,Tgas,FlashFactor,h_oil_meas,PAHH,PALL,Oil_LAHH,Oil_LALL,Orientation,Ts,Tag); 
    
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                                
    sys = mdlOutputs(t,x,u,rho_oil,h_oil_meas);

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case {1,4,9}, % unused flags                  
    sys = []; % do nothing

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    error(['unhandled flag = ',num2str(flag)]);
end

%
%=======================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=======================================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes(R_sep,L_sep,n,ho0,p0,Orientation)

sizes = simsizes;
sizes.NumContStates  = 0;
sizes.NumDiscStates  = 4;
sizes.NumOutputs     = 3;
sizes.NumInputs      = 3;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

    
V_sep=R_sep^2*pi*L_sep+(pi/4)*(4/3)*R_sep^3;    
%caluclate initial oil volume
if Orientation==1 %Horizontal orientation
    Vo0 = OilVol(ho0,R_sep,L_sep,n);
else %Vertical orientation
    Vo0=(pi/4)*(2/3)*R_sep^3+R_sep^2*pi*ho0;
end
Vg0 = V_sep - Vo0;

if ((Vg0 < 0) | (Vg0 > V_sep))
   errordlg('Something wrong with initial gas volume');
end 
x0  = [ho0 Vg0 p0*1e5 Vo0]; % initialize states
str = [];
ts  = [-1 0]; 

% end mdlInitializeSizes
function sys = mdlUpdate(t,x,u,R,rho_gasS,MW_gas,rho_oilS,rho_oil,R_sep,L_sep,n,Tgas,FlashFactor,h_oil_meas,PAHH,PALL,Oil_LAHH,Oil_LALL,Orientation,Ts,Tag)
persistent errorStatusOil errorStatusPressure
if ~isfield(errorStatusOil,Tag)
    errorStatusOil=setfield(errorStatusOil,Tag,0)
    errorStatusPressure=setfield(errorStatusPressure,Tag,0)
end
% Convert inputs
n_gas_in=u(1)*rho_oilS*FlashFactor/100/MW_gas; %Sm3/s oil -> kg/s oil -> kg/s gas (flashed off) ->mol/s gas
n_gas_out = u(2)*rho_gasS/MW_gas; % Sm3/s -> mol/s
Q_oil_in = u(1)*rho_oilS*(100-FlashFactor)/100/rho_oil; % Sm3/s -> m3/s
Q_oil_out = u(3)*rho_oilS/rho_oil; %Sm3/s -> m3/s

% options = optimset('LargeScale','off');
options = optimset('LargeScale','off','Display','off');

% oil level
DeltaVo = (Q_oil_in-Q_oil_out)*Ts;
if Orientation==1 %Horizontal orientation
    sys(1) = fzero(@hliq,[x(1)],options,R_sep,L_sep,n,x(4)+DeltaVo);
else %Vertical orientation
    sys(1)=x(1)+DeltaVo/(pi*R_sep^2);
end

if isreal(sys(1)) == 0
    disp('Oil volume: Complex solution');
end

if sys(1) <= Oil_LALL & ~getfield(errorStatusOil,Tag),
    errordlg('Trip on oil level (LALL)');
    errorStatusOil=setfield(errorStatusOil,Tag,1);
elseif sys(1) >= Oil_LAHH & ~getfield(errorStatusOil,Tag),
    errordlg('Trip on oil level (LAHH)');
    errorStatusOil=setfield(errorStatusOil,Tag,1);
elseif getfield(errorStatusOil,Tag) & (sys(1) < Oil_LAHH) & (sys(1) > Oil_LALL)
    errorStatusOil=setfield(errorStatusOil,Tag,0);
end

% Gas volume
sys(4) = x(4)+DeltaVo;
V_sep=R_sep^2*pi*L_sep+(pi/4)*(4/3)*R_sep^3;
sys(2) = V_sep-(x(4)+DeltaVo);
if ((sys(2) < 0) | (sys(2) > V_sep))
    errordlg('Something wrong with calculated gas volume');
end
% Gas pressure
sys(3) = x(3) + ((n_gas_in-n_gas_out)*R*Tgas*Ts)/x(2)...
    -(Ts*x(3)/x(2))*(-Q_oil_in+Q_oil_out);
if sys(3) <= PALL*1e5 & ~getfield(errorStatusPressure,Tag),
    errordlg('Trip on gas pressure (PALL)');
    errorStatusPressure=setfield(errorStatusPressure,Tag,1);
elseif sys(3) >= PAHH*1e5 & ~getfield(errorStatusPressure,Tag),
    errordlg('Trip on gas pressure (PAHH)');
    errorStatusPressure=setfield(errorStatusPressure,Tag,1);
elseif getfield(errorStatusPressure,Tag) & (sys(3) < PAHH*1e5) & (sys(3) > PALL*1e5)
    errorStatusPressure=setfield(errorStatusPressure,Tag,0);
end
% end mdlUpdate
%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u,rho_oil,h_oil_meas)
sys(1) = x(3)/(1e5); % Pa -> bara/bar
sys(2) = x(1); % Oil level above bottom
sys(3) = (x(3) + rho_oil*9.81*(x(1)-h_oil_meas))/(1e5); %SIMPLIFICATION, IGNORING WATER CONTENT
% end mdlOutputs


function OilVol = OilVol(ho,R,L,n) 
OilVol = cyl(ho,R,L)+ 2*one_ellip(ho,R,n);


function cyl=cyl(h,R,L)
if h<0
    h=0;
elseif h>2*R
    h=2*R;
end;
root = 2*R*h-h^2;
cyl  = L*((h-R)*sqrt(root)+R^2*(asin((h-R)/R)+pi/2));

function one_ellip=one_ellip(h,R,n)
one_ellip = (pi/(2*n*n))*(R*h*h-(h*h*h)/3);

function f = hliq(ho,R,L,n,Vo)
f = cyl(ho,R,L)+2*one_ellip(ho,R,n)-Vo;
f = (f/Vo);

