function [sys,x0,str,ts] = pwflashdrum(t,x,u,flag,R,rho_gasS,MW_gas,rho_oilS,...
    rho_oil,rho_waterS,rho_water,R_sep,L_sep,n,hweir_min,hweir_max,Lweir,Tgas...
    ,h_oil_meas,hw0,ho0,p0,PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,Tag,wcc0,wcd0)


% 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
% hweir               Weir height [m]
% Lweir               Length to weir [m]
% 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]
% W_LAHH              High level trip, water [m]
% W_LALL              Low level trip, water [m]
%
% 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. Oil rate out [Sm3/s]
% 6. Water rate out [Sm3/s]
%
% Outputs:
% 1. Gas pressure [bara]
% 2. Oil level [m]
% 3. Water level [m]
% 4. Pressure @ oil output [bara]
% 5. Pressure @ water output [bara]


switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts] = mdlInitializeSizes(R_sep,L_sep,n,hweir_max,Lweir,hw0,ho0,p0,wcc0,wcd0);
  %%%%%%%%%%%%%%%%%%%%%%%%%
  % Computing 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,hweir_min,hweir_max,Lweir,Tgas,...
    PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,Tag); 
    
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                                
    sys = mdlOutputs(t,x,u,rho_oil,rho_water,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,...
    hweir_max,Lweir,hw0,ho0,p0,wcc0,wcd0)

sizes = simsizes;
sizes.NumContStates  = 8;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 7;
sizes.NumInputs      = 8; 
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

%calculate initial volum on clean side
Vclean = WaterVol(hw0,R_sep,Lweir,L_sep,n,hweir_max);    
Vwclean=Vclean*wcc0;
Voclean=Vclean*(1-wcc0);

%caluclate initial oil volume
Vdirty = OilVol(ho0,R_sep,Lweir,L_sep,n,hweir_max,Vclean);
Vwdirty=Vdirty*wcd0;
Vodirty=Vdirty*(1-wcd0);
V_sep=R_sep^2*pi*L_sep+(pi/4)*(4/3)*R_sep^3;
Vg0 = V_sep - Vclean - Vdirty;

if ((Vg0 < 0) | (Vg0 > V_sep))
   errordlg('Something wrong with initial gas volume');
end 
x0  = [hw0 ho0 Vg0 p0*1e5 Vwclean Vwdirty Voclean Vodirty]; % 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,hweir_min,hweir_max,Lweir,Tgas,...
    PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,Tag)
% Convert inputs
persistent errorStatusOil errorStatusWater errorStatusPressure
if ~isfield(errorStatusOil,Tag)
    errorStatusOil=setfield(errorStatusOil,Tag,0);
    errorStatusWater=setfield(errorStatusWater,Tag,0);
    errorStatusPressure=setfield(errorStatusPressure,Tag,0);
end
   
% -
% --------------------------------------------------------------

% Inputs:
% 1. Total gas rate in [Sm3/s]
% 2. Liquid rate in dirty inflow [Sm3/s]
% 3. Water cut in dirty inflow [0-1]
% 4. Liquid rate in clean inflow [Sm3/s]
% 5. Water cut in dirty inflow [0-1]
% 6. Gas rate out [Sm3/s]
% 7. Dirty rate out [Sm3/s]
% 8. Clean rate out [Sm3/s]

% Water cuts:
wcd=x(6)/(x(6)+x(8));
wcc=x(5)/(x(5)+x(7));

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_dirty = (u(2)*(1-u(3)))*rho_oilS/rho_oil; % Sm3/s -> m3/s
Q_oil_in_clean = (u(4)*(1-u(5)))*rho_oilS/rho_oil; % Sm3/s -> m3/s
Q_oil_out_dirty = u(7)*rho_oilS/rho_oil*(1-wcd); % Sm3/s -> m3/s
Q_oil_out_clean = u(8)*rho_oilS/rho_oil*(1-wcc); % Sm3/s -> m3/s 

Q_water_in_dirty = u(2)*u(3)*rho_waterS/rho_water; % Sm3/s -> m3/s
Q_water_in_clean = u(4)*u(5)*rho_waterS/rho_water; % Sm3/s -> m3/s
Q_water_out_dirty = u(7)*rho_waterS/rho_water*wcd;  % Sm3/s -> m3/s
Q_water_out_clean = u(8)*rho_waterS/rho_water;  % Sm3/s -> m3/s

VWeir_min = cyl(hweir_min,R_sep,Lweir) + one_ellip(hweir_min,R_sep,n);
Q_overflow=max((x(1)-hweir_min)/(hweir_max-hweir_min)*(x(5)+x(7)-VWeir_min),0)*0.0001;
% assuming perfect separation on clean side
if Q_overflow>0
    wc_overflow=max((x(5)-VWeir_min)/(x(5)+x(7)),0);
    Q_overflow_oil=Q_overflow*(1-wc_overflow);
    Q_overflow_water=Q_overflow*(wc_overflow);
else
    Q_overflow_oil=0;
    Q_overflow_water=0;
end

% clean side
dVocdt=Q_oil_in_clean-Q_overflow_oil-Q_oil_out_clean;
dVwcdt=Q_water_in_clean-Q_overflow_water-Q_water_out_clean;
dVcdt = dVocdt+dVwcdt;
dVcdhw=2*Lweir*sqrt(2*R_sep*x(1)-x(1)^2)+pi/2/n^2*(2*R_sep*x(1)-x(1)^2);
dhwdt=inv(dVcdhw)*dVcdt;
if x(1) <= W_LALL & ~getfield(errorStatusWater,Tag),
    errordlg(['Trip on water level (LALL) in ' Tag] );
    errorStatusWater=setfield(errorStatusWater,Tag,1);
elseif x(1) >= W_LAHH & ~getfield(errorStatusWater,Tag),
    errordlg(['Trip on water level (LAHH) in ' Tag]);
    errorStatusWater=setfield(errorStatusWater,Tag,1);
elseif getfield(errorStatusWater,Tag) & (x(1) < W_LAHH) & (x(1) > W_LALL)
    errorStatusWater=setfield(errorStatusWater,Tag,0);
end
% dirty side
dVoddt=Q_oil_in_dirty+Q_overflow_oil-Q_oil_out_dirty;
dVwddt=Q_water_in_dirty+Q_overflow_water-Q_water_out_dirty;
dVddt = dVoddt+dVwddt;
dVddhw=2*(L_sep-Lweir)*sqrt(2*R_sep*x(2)-x(2)^2)+pi/2/n^2*(2*R_sep*x(2)-x(2)^2);
dhodt=inv(dVddhw)*(dVddt);

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

% Gas volume
dVgdt=-(dVcdt+dVddt);
% Gas pressure
dngdt=n_gas_in-n_gas_out;
dPdt=R*Tgas/x(3)*dngdt-x(4)/x(3)*dVgdt;
if x(4) <= PALL*1e5 & ~getfield(errorStatusPressure,Tag),
    errordlg(['Trip on gas pressure (PALL) in ' Tag]);
    errorStatusPressure=setfield(errorStatusPressure,Tag,1);
elseif x(4) >= PAHH*1e5 & ~getfield(errorStatusPressure,Tag),
    errordlg(['Trip on gas pressure (PAHH) in ' Tag]);
    errorStatusPressure=setfield(errorStatusPressure,Tag,1);
elseif getfield(errorStatusPressure,Tag) & (x(4) < PAHH*1e5) & (x(4) > PALL*1e5)
    errorStatusPressure=setfield(errorStatusPressure,Tag,0);
end
sys=[dhwdt dhodt dVgdt dPdt dVwcdt dVwddt dVocdt dVoddt];
% end mdlDerivatives
%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u,rho_oil,rho_water,h_oil_meas)
wcd=x(6)/(x(6)+x(8));
wcc=x(5)/(x(5)+x(7));

sys(1) = x(4)/(1e5); % Pa -> bara/bar
sys(2) = x(2); % level dirty side
sys(3) = x(1); % level clean side
sys(4) = (x(4) + rho_oil*9.81*(x(2)-h_oil_meas))/(1e5); %pressure @ dirty outlet
sys(5) = (x(4) + rho_water*9.81*x(1))/(1e5); %pressure @ clean outlet
sys(6) = wcd;
sys(7) = wcc;

% end mdlOutputs

function WaterVol = WaterVol(hw,R,Lweir,L,n,hweir)
if hw <= hweir % water below weir
    WaterVol = cyl(hw,R,Lweir) + one_ellip(hw,R,n);
else % water above weir -> This is not valid
    errordlg('Water above weir');
    WaterVol = cyl(hw,R,Lweir) + one_ellip(hw,R,n);
end


function OilVol = OilVol(ho,R,Lweir,L,n,hweir,Vw0) % Excluding oil on water side
if ho <= hweir % oil below weir
    OilVol = cyl(ho,R,L-Lweir) + one_ellip(ho,R,n);
else 
    OilVol = cyl(ho,R,L) + 2*one_ellip(ho,R,n) ... % as if only oil is in the separator
            - Vw0; % Subtract water
%    OilVol = cyl(ho,R,L)-cyl(hweir,R,Lweir) + 2*one_ellip(ho,R,n) - one_ellip(hweir,R,n); % Disregard water side
end

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);



