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

% For simplicity, this model assumes that oil is always filling up the
% water side (up to the weir height)


% 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. Gas rate in [Sm3/s]
% 2. Oil rate in [Sm3/s]
% 3. Water rate in [Sm3/s]
% 4. Gas rate out [Sm3/s]
% 5. Liquid 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);
  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,                                                
    sys = mdlDerivatives(t,x,u,R,rho_gasS,MW_gas,rho_oilS,rho_oil,rho_waterS,rho_water,R_sep,L_sep,n,Tgas,h_oil_meas,PAHH,PALL,Oil_LAHH,Oil_LALL,Orientation,Tag); 
    
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                                
    sys = mdlOutputs(t,x,u,rho_oil,h_oil_meas);

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case {2,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  = 4;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 3;
sizes.NumInputs      = 5; 
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

    
    
%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
V_sep=R_sep^2*pi*L_sep+(pi/4)*(4/3)*R_sep^3;
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  = [0 0]; 

% end mdlInitializeSizes
function sys = mdlDerivatives(t,x,u,R,rho_gasS,MW_gas,rho_oilS,rho_oil,rho_waterS,rho_water,R_sep,L_sep,n,Tgas,h_oil_meas,PAHH,PALL,Oil_LAHH,Oil_LALL,Orientation,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_gasS/MW_gas; % Sm3/s -> mol/s
n_gas_out = u(4)*rho_gasS/MW_gas; % Sm3/s -> mol/s
Q_oil_in = u(2)*rho_oilS/rho_oil; % Sm3/s -> m3/s
Q_liq_out = u(5);
Q_water_in = u(3)*rho_waterS/rho_water; % Sm3/s -> m3/s
Q_liq_in=Q_oil_in+Q_water_in;

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

% oil level
dVldt = (Q_liq_in-Q_liq_out);
if Orientation==1 %Horizontal orientation
    dVldhl=2*L_sep*sqrt(2*R_sep*x(1)-x(1)^2)+pi/n^2*(2*R_sep*x(1)-x(1)^2);
    dhodt=inv(dVldhl)*(dVldt);
else %Vertical orientation
    dhodt=dVldt/(pi*R_sep^2);
end

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


% Gas volume
dVgdt=-dVldt;
V_sep=R_sep^2*pi*L_sep+(pi/4)*(4/3)*R_sep^3;
% Gas pressure
dngdt=n_gas_in-n_gas_out;
dPdt=R*Tgas/x(2)*dngdt-x(3)/x(2)*dVgdt;
if x(3) <= PALL*1e5 & ~getfield(errorStatusPressure,Tag),
    errordlg('Trip on gas pressure (PALL)');
    errorStatusPressure=setfield(errorStatusPressure,Tag,1);
elseif x(3) >= PAHH*1e5 & ~getfield(errorStatusPressure,Tag),
    errordlg('Trip on gas pressure (PAHH)');
    errorStatusPressure=setfield(errorStatusPressure,Tag,1);
elseif getfield(errorStatusPressure,Tag) & (x(3) < PAHH*1e5) & (x(3) > PALL*1e5)
    errorStatusPressure=setfield(errorStatusPressure,Tag,0);
end
sys=[dhodt dVgdt dPdt dVldt];
% 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);

