function [sys,x0,str,ts] = sSeparator(t,x,u,flag,R,rho_gasS,MW_gas, ...
    rho_oilS,rho_oil,rho_waterS,rho_water,R_sep,L_sep,n,hweir,Lweir,Tgas, ...
    h_oil_meas,hw0,ho0,p0,PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,SepTag, ...
    wc0,delta_t,N_drop,q_mean_0,tau_0)

% This s-function present a 3 phase separator including a weir separating
% the two liquid outlets. The model is dynamic.
%
% For simplicity, this model assumes that oil is always filling up the
% water side (up to the weir height)

global wc

% 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]
% SepTag              Tag ('String') Separator name/number
% wc0                 Initial water cut, fraction between 0 and 1
% delta_t             Sampling time in Q
% N_drop              Efficiency of the separator, that is how much oil
% goes from the water phase to the oil phase.
% 
% 
% 
% Inputs:
% All of the in and out flows are physical ones, that is water in and out
% includes the disperged oil.
% 1. Gas rate in [Sm3/s]
% 2. Oil rate in [Sm3/s]
% 
% The oil fractions are fractions of disperged oil in the water phase
% 3. x_oilInDispSep_1, fraction of oil in chamber 1
% 4. x_oilInDispSep_2, fraction of oil in chamber 2
% 5. x_oilInDispSep_3, fraction of oil in chamber 3
% 6. x_oilInDispSep_4, fraction of oil in chamber 4
% 7. x_oilInDispSep_5, fraction of oil in chamber 5

% 8. Water rate in [Sm3/s] 
% 9 Gas rate out [Sm3/s]
% 10. Oil rate out [Sm3/s] 
% 11. 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]This will also be the pressure going
% into the hydrocyclone.
% 6. x_oilOutDispSep_1, Fraction of oil out with the water in chamber 1
% 7. x_oilOutDispSep_2, Fraction of oil out with the water in chamber 2
% 8. x_oilOutDispSep_3, Fraction of oil out with the water in chamber 3
% 9. x_oilOutDispSep_4, Fraction of oil out with the water in chamber 4
% 10. x_oilOutDispSep_5, Fraction of oil out with the water in chamber 5
% 11. x_oilOutDispSep_tot, total oil fraction out with the water
% 12. tau, the retention time of the disperged oil in the separator
% 13- tau_cstr. The retentiontime for the separator


switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts] = mdlInitializeSizes(R_sep,L_sep,n,hweir,Lweir,hw0,ho0, ...
         p0,wc0,q_mean_0);
  %%%%%%%%%%%%%%%%%%%%%%%%%
  % 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,Lweir,Tgas,h_oil_meas, ...
        PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,SepTag,wc0,hw0,delta_t,N_drop,tau_0);
    
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                                
    sys = mdlOutputs(t,x,u,flag,R,rho_gasS,MW_gas,rho_oilS,...
    rho_oil,rho_waterS,rho_water,R_sep,L_sep,n,hweir,Lweir,Tgas,h_oil_meas, ...
    hw0,ho0,p0,PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,SepTag,wc0,delta_t, ...
    N_drop,tau_0);

  %%%%%%%%%%%%%
  % 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,Lweir,hw0, ...
     ho0,p0,wc0,q_mean_0) 
     
global wc

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

sys = simsizes(sizes);
%calculate initial water volume

if hw0 <= hweir
    Vw0 = WaterVol(R_sep,Lweir,L_sep,n,hweir,hw0);
else
    errordlg('Water level above weir plate');
end
    
wc = wc0;
    
%caluclate initial oil volume
VwWeir = cyl(hweir,R_sep,Lweir) + one_ellip(hweir,R_sep,n);
Vo0 = OilVol(ho0,R_sep,Lweir,L_sep,n,hweir,Vw0);
V_sep=R_sep^2*pi*L_sep+(pi/4)*(4/3)*R_sep^3;
Vg0 = V_sep - Vo0 - Vw0;

if ((Vg0 < 0) | (Vg0 > V_sep))
   errordlg('Something wrong with initial gas volume');
end 
x0  = [hw0 ho0 Vg0 p0*1e5 Vw0 Vo0 (Vw0/q_mean_0)+10]; % 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,Lweir,Tgas,h_oil_meas,PAHH, ...
    PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,SepTag,wc0,hw0,delta_t,N_drop,tau_0)
% Convert inputs
persistent errorStatusOil errorStatusWater errorStatusPressure
global wc

if ~isfield(errorStatusOil,SepTag)
    errorStatusOil=setfield(errorStatusOil,SepTag,0);
    errorStatusWater=setfield(errorStatusWater,SepTag,0);
    errorStatusPressure=setfield(errorStatusPressure,SepTag,0);
end

% calculate water volume if water level is equal to weir height:
VwWeir = cyl(hweir,R_sep,Lweir) + one_ellip(hweir,R_sep,n);
% calculate oil volume if oil level is equal to weir height:
VoWeir = cyl(hweir,R_sep,L_sep-Lweir) + one_ellip(hweir,R_sep,n);

% --------------------------------------------------------------
% - modification @ 14.03.06 - ER
%
% - calculate total liquid volume
if (x(1)>x(2))  % h_water > h_oil
    hl = x(1);
else            % h_oil >= h_water
    hl = x(2);
end
Vliq = cyl(hl,R_sep,L_sep) + 2*one_ellip(hl,R_sep,n);

% - calculate water cut
if (wc0 < 0 || wc0 > 1)
    errordlg(['Input water cut out of range in ' SepTag] );
end
wc = wc0;

if ( x(1) > hweir && x(5) < VwWeir)
    disp ('no')
end

if ( x(5) > VwWeir )
    wc = (x(5) - VwWeir + wc0*x(6))/(Vliq-VwWeir);
end

if (wc < 0 || wc > 1)
    errordlg(['Calculated water cut out of range in ' SepTag] );
end

n_gas_in = u(1)*rho_gasS/MW_gas; % Sm3/s -> mol/s
n_gas_out = u(9)*rho_gasS/MW_gas; % Sm3/s -> mol/s

Q_oil_in = u(2)*rho_oilS/rho_oil; % Sm3/s -> m3/s

Q_water_in = u(8)*rho_waterS/rho_water; % Sm3/s -> m3/s

%x_oilInSep_0, This is the fraction which will be building up, and go
% over to the oil phase at the end.



% Added the oil separated out from the water phase
% modified with water cut - 14.03.06 - ER
if t==0
Q_oilOutWoil_tot=0;
end
if t==0
Q_water_out=Q_water_in;

else
[Q_water_out,Q_oilOutWoil_tot,x_oilOutDispSep_1,x_oilOutDispSep_2, ...
     x_oilOutDispSep_3,x_oilOutDispSep_4,x_oilOutDispSep_5, ...
      x_oilOutDispSep_tot,tau]=WaterDisp(x,u,rho_oilS,rho_oil,rho_waterS, ...
       rho_water,R_sep,Lweir,L_sep,n,hweir,x(1),delta_t,t,N_drop,tau_0);
%Q_water_out = u(11)*rho_waterS/rho_water + u(2)*rho_oilS/rho_oil*wc-Q_oilOutWoil_tot; % Sm3/s -> m3/s
% Subtracting the oil separated out
end

Q_oil_out = u(10)*rho_oilS/rho_oil*(1-wc);%+(0.2*Q_water_in-Q_oilOutSep_tot); % Sm3/s -> m3/s
% modified with water cut - 14.03.06 - ER


% options = optimset('LargeScale','off');
%options = optimset('Display','off');
q_mean= (Q_water_in+Q_water_out)/2;
dqmeandt=u(12);
% water volume
dVwdt = (Q_water_in-Q_water_out);
dVwdhw=2*Lweir*sqrt(2*R_sep*x(1)-x(1)^2)+pi/2/n^2*(2*R_sep*x(1)-x(1)^2);
dTwdt=(dVwdt*q_mean-x(5)*dqmeandt)/q_mean^2;
dhwdt=inv(dVwdhw)*dVwdt;
if x(1) <= W_LALL & ~getfield(errorStatusWater,SepTag),
    errordlg(['Trip on water level (LALL) in ' SepTag] );
    errorStatusWater=setfield(errorStatusWater,SepTag,1);
elseif x(1) >= W_LAHH & ~getfield(errorStatusWater,SepTag),
    errordlg(['Trip on water level (LAHH) in ' SepTag]);
    errorStatusWater=setfield(errorStatusWater,SepTag,1);
elseif getfield(errorStatusWater,SepTag) & (x(1) < W_LAHH) & (x(1) > W_LALL)
    errorStatusWater=setfield(errorStatusWater,SepTag,0);
end
% oil level
dVodt = (Q_oil_in-Q_oil_out+Q_oilOutWoil_tot);
dVldhw=2*L_sep*sqrt(2*R_sep*x(2)-x(2)^2)+pi/n^2*(2*R_sep*x(2)-x(2)^2);
dhodt=inv(dVldhw)*(dVwdt+dVodt);

if x(2) <= Oil_LALL & ~getfield(errorStatusOil,SepTag),
    errordlg(['Trip on oil level (LALL) in ' SepTag]);
    errorStatusOil=setfield(errorStatusOil,SepTag,1);
elseif x(2) >= Oil_LAHH & ~getfield(errorStatusOil,SepTag),
    errordlg(['Trip on oil level (LAHH) in' SepTag]);
    errorStatusOil=setfield(errorStatusOil,SepTag,1);
elseif getfield(errorStatusOil,SepTag) & (x(2) < Oil_LAHH) & (x(2) > Oil_LALL)
    errorStatusOil=setfield(errorStatusOil,SepTag,0);
end
% Gas volume
dVgdt=-(dVwdt+dVodt);
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(3)*dngdt-x(4)/x(3)*dVgdt;
if x(4) <= PALL*1e5 & ~getfield(errorStatusPressure,SepTag),
    errordlg(['Trip on gas pressure (PALL) in ' SepTag]);
    errorStatusPressure=setfield(errorStatusPressure,SepTag,1);
elseif x(4) >= PAHH*1e5 & ~getfield(errorStatusPressure,SepTag),
    errordlg(['Trip on gas pressure (PAHH) in ' SepTag]);
    errorStatusPressure=setfield(errorStatusPressure,SepTag,1);
elseif getfield(errorStatusPressure,SepTag) & (x(4) < PAHH*1e5) & (x(4) > PALL*1e5)
    errorStatusPressure=setfield(errorStatusPressure,SepTag,0);
end
sys=[dhwdt dhodt dVgdt dPdt dVwdt dVodt dTwdt];
% end mdlDerivatives
%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u,flag,R,rho_gasS,MW_gas,rho_oilS,...
    rho_oil,rho_waterS,rho_water,R_sep,L_sep,n,hweir,Lweir,Tgas,h_oil_meas, ...
    hw0,ho0,p0,PAHH,PALL,Oil_LAHH,Oil_LALL,W_LAHH,W_LALL,SepTag,wc0,delta_t, ...
    N_drop,tau_0)
    
%% Including a disperged oil in water phase, with a separation efficiency
%% depending on the retention time and coallecence following second order
%% dynamics



% Outputs:
% 1. Gas pressure [bara]
% 2. Oil level [m]
% 3. Water level [m]
% 4. Pressure @ oil output [bara]
% 5. Pressure @ water output [bara]
% 6. x_oilOutDispSep_1, Fraction of oil out with the water in chamber 1
% 7. x_oilOutDispSep_2, Fraction of oil out with the water in chamber 2
% 8. x_oilOutDispSep_3, Fraction of oil out with the water in chamber 3
% 9. x_oilOutDispSep_4, Fraction of oil out with the water in chamber 4
% 10. x_oilOutDispSep_5, Fraction of oil out with the water in chamber 5
% 11. x_oilOutDispSep_tot, total oil fraction out with the water
% 12. tau, the retention time of the disperged oil in the separator
 
[Q_water_out,Q_oilOutWoil_tot,x_oilOutDispSep_1,x_oilOutDispSep_2, ... 
    x_oilOutDispSep_3,x_oilOutDispSep_4,x_oilOutDispSep_5, ... 
    x_oilOutDispSep_tot,tau]=WaterDisp(x,u,rho_oilS,rho_oil,rho_waterS, ... 
    rho_water,R_sep,Lweir,L_sep,n,hweir,x(1),delta_t,t,N_drop,tau_0);

sys(1) = x(4)/(1e5); % Pa -> bara/bar
sys(2) = x(2); % Oil level above bottom
sys(3) = x(1);
sys(4) = (x(4) + rho_oil*9.81*(x(2)-h_oil_meas))/(1e5);
sys(5) = (x(4) + rho_oil*9.81*(x(2)-x(1)) + rho_water*9.81*x(1))/(1e5);
sys(6) = x_oilOutDispSep_1;
sys(7) = x_oilOutDispSep_2;
sys(8) = x_oilOutDispSep_3;
sys(9) = x_oilOutDispSep_4;
sys(10) = x_oilOutDispSep_5;
sys(11) = x_oilOutDispSep_tot;
sys(12) = tau;
sys(13) = x(5)/u(8);%x(7);

% try
%     sys(6) = wc;
% catch
%     sys(6) = wc0;
% end

% end mdlOutputs

function [Q_water_out,Q_oilOutWoil_tot,x_oilOutDispSep_1,x_oilOutDispSep_2, ...
     x_oilOutDispSep_3,x_oilOutDispSep_4,x_oilOutDispSep_5, ... 
     x_oilOutDispSep_tot,tau]=WaterDisp(x,u,rho_oilS,rho_oil,rho_waterS, ... 
     rho_water,R_sep,Lweir,L_sep,n,hweir,hw0,delta_t,t,N_drop,tau_0)
% Inputs:
% All of the in and out flows are physical ones, that is water in and out
% includes the disperged oil.
% 1. Gas rate in [Sm3/s]
% 2. Oil rate in [Sm3/s]
% 
% The oil fractions are fractions of disperged oil in the water phase
% 4. x_oilInSep_1, fraction of oil in chamber 1
% 5. x_oilInSep_2, fraction of oil in chamber 2
% 6. x_oilInSep_3, fraction of oil in chamber 3
% 7. x_oilInSep_4, fraction of oil in chamber 4
% The fraction of oil in chamber 5 will given by the total oil mass balance
% 
% 8. Water rate in [Sm3/s] 
% 9 Gas rate out [Sm3/s]
% 10. Oil rate out [Sm3/s] 
% 11. Water rate out [Sm3/s]
%  
global wc

%Defining the inputs
x_oilInDispSep_1 = u(3);
x_oilInDispSep_2 = u(4);
x_oilInDispSep_3 = u(5);
x_oilInDispSep_4 = u(6);
x_oilInDispSep_5 = u(7);

Q_water_in = u(8)*rho_waterS/rho_water; % Sm3/s -> m3/s
Q_oil_in = u(2)*rho_oilS/rho_oil; % Sm3/s -> m3/s

% Calculating the fraction of the remaining chamber
x_oilInDispSep = x_oilInDispSep_1 + x_oilInDispSep_2 + x_oilInDispSep_3 + ...
     x_oilInDispSep_4+x_oilInDispSep_5; 
Q_oilInDispSep=Q_water_in*x_oilInDispSep;

x_oilInDispSep_vec = [x_oilInDispSep_1 x_oilInDispSep_2 x_oilInDispSep_3 ... 
    x_oilInDispSep_4 x_oilInDispSep_5];

%Finding the different compartment flows
QoilDispSep_in1=x_oilInDispSep_1*Q_water_in; %m3/s
QoilDispSep_in2=x_oilInDispSep_2*Q_water_in; %m3/s
QoilDispSep_in3=x_oilInDispSep_3*Q_water_in; %m3/s
QoilDispSep_in4=x_oilInDispSep_4*Q_water_in;%m3/s
QoilDispSep_in5=x_oilInDispSep_5*Q_water_in; %m3/s

Q_chambers = [QoilDispSep_in1 QoilDispSep_in2 QoilDispSep_in3 ...
     QoilDispSep_in4 QoilDispSep_in5];

% Introducing 2nd order characteristics for the plug flow going through
% Given by (1/CA)=(1/CA0)+kt. Assuming CA = x_oilInSep and t=tau and that k
% is given by assuming complete separation after 4 minutes, and that CA at
% thi point equals 99.9 % of CA0.

k1 = ((1/(0.001*x_oilInDispSep_1))-(1/x_oilInDispSep_1))/(4*60);
k2 = ((1/(0.001*x_oilInDispSep_2))-(1/x_oilInDispSep_2))/(4*60);
k3 = ((1/(0.001*x_oilInDispSep_3))-(1/x_oilInDispSep_3))/(4*60);
k4 = ((1/(0.001*x_oilInDispSep_4))-(1/x_oilInDispSep_4))/(4*60);
k5 = ((1/(0.001*x_oilInDispSep_5))-(1/x_oilInDispSep_5))/(4*60);

k_vec = [k1 k2 k3 k4 k5];

persistent Qvec

if t == 0
    Qvec = Q_water_in;
else
    Qvec = [Qvec Q_water_in];
end

% What is denoted a/A will be going upwards/be separated out

Qoil_a = 0;

% What is denoted c/C on the other hand will be going
% downward/exiting still disperged in the water at the outlet.

Qoil_c = 0;

Q_oilOutDispSep_1 = 0;
Q_oilOutDispSep_2 = 0;
Q_oilOutDispSep_3 = 0;
Q_oilOutDispSep_4 = 0;
Q_oilOutDispSep_5 = 0;

Q_oilOutWOil_1 = 0;
Q_oilOutWOil_2 = 0;
Q_oilOutWOil_3 = 0;
Q_oilOutWOil_4 = 0;
Q_oilOutWOil_5 = 0;

for i=1:length(Q_chambers)
    
if t<=250  
    
    %Tanken er at scriptet bør kjøres en tid først (>= en oppholdstid
    %slik at Q vektoren får tid og bygge seg opp) Lagre verdier bakover i 
    %tid vha persistent. See the residence time function further down in
    %the script
    
    % First the minimun efficiency based on drop size is estimated
    Qoil_a(i)=N_drop(i)*Q_chambers(i);
    Qoil_c(i)=(1-N_drop(i))*Q_chambers(i);
    
    tau=tau_0; % An initial retentiontime is provided for the initialisation 
            %Adding an initial residence time dependant on the script 
            %running, put out as parameter in mask?

            %Succedingly the separation provided from the retentiontime the
            %fluid spends in the separatior comes in addition.
            
        x_oilOutDispSep(i) = x_oilInDispSep_vec(i)*(1-N_drop(i))/ ...
             (1+k_vec(i)*tau*x_oilInDispSep_vec(i)*(1-N_drop(i)));  
        %Denne er i vann
        
       
        
         QoilA_tau(i)=(x_oilInDispSep_vec(i)-x_oilOutDispSep(i))*Q_chambers(i);
         QoilC_tau(i)=x_oilOutDispSep(i)*Q_chambers(i);
       
         % Placing the drops in their respective compartments
       if i==1;
            Q_oilOutWOil_1 = QoilA_tau(i);
            Q_oilOutDispSep_1 = QoilC_tau(i);
        elseif i==2
            Q_oilOutWOil_2 = QoilA_tau(i);
            Q_oilOutDispSep_2 = QoilC_tau(i);
        elseif i==3
            Q_oilOutWOil_3 = QoilA_tau(i);
            Q_oilOutDispSep_3 = QoilC_tau(i);
        elseif i==4
            Q_oilOutWOil_4 = QoilA_tau(i);
            Q_oilOutDispSep_4 = QoilC_tau(i);
        elseif i==5
            Q_oilOutWOil_5 = QoilA_tau(i);
            Q_oilOutDispSep_5 = QoilC_tau(i);
        end

else %Etter en startkjøring begynner å ta hensyn til oppholdstidens endring
    %med tiden.
  
tau = residencetime(R_sep,Lweir,L_sep,n,hweir,hw0,Qvec,delta_t);  
  
    Qoil_a(i)=N_drop(i)*Q_chambers(i); % Efferktivitet på dråpestr basis 
    % Q_a er det som går opp til emulsjonslaget
    Qoil_c(i)=(1-N_drop(i))*Q_chambers(i);
    
    if tau<60
        NA_tau=0; %0.2;
        NC_tau=1; %0.5;
        
        %Det du egentlig sier her er at ved en for liten oppholdstid vil
        %ingenting gå opp... Men du må bestemme deg for om dette skal
        %gjelde absolutt eller om de fremdeles vil ha effekten gitt av
        %dråpestr... 
        
       QoilA_tau(i)=0;
        QoilC_tau(i)=Qoil_a(i)+Qoil_c(i);
        
    elseif tau<240
        % NA_tau er her fraksjonen som går opp (separeres fra) Husk å
        % definer dette klart. Er egentlig CA, og CA = CA0/(1+k*t(au)*CA0)
        
        x_oilOutDispSep(i) = x_oilInDispSep_vec(i)*(1-N_drop(i))/ ...
             (1+k_vec(i)*tau*x_oilInDispSep_vec(i)*(1-N_drop(i)));   
        %Vannfraksjon
        
        QoilA_tau(i)=(x_oilInDispSep_vec(i)-x_oilOutDispSep(i))*Q_chambers(i);
         QoilC_tau(i)=x_oilOutDispSep(i)*Q_chambers(i);
       
    else
        NA_tau=1; 
        NC_tau=0;
        
        QoilA_tau(i)=Qoil_a(i)+Qoil_c(i);
        QoilC_tau(i)=0;
    end

     % Placing the drops in their respective compartments

        if i==1;
            Q_oilOutWOil_1 = QoilA_tau(i);
            Q_oilOutDispSep_1 = QoilC_tau(i);
        elseif i==2
            Q_oilOutWOil_2 = QoilA_tau(i);
            Q_oilOutDispSep_2 = QoilC_tau(i);
        elseif i==3
            Q_oilOutWOil_3 = QoilA_tau(i);
            Q_oilOutDispSep_3 = QoilC_tau(i);
        elseif i==4
            Q_oilOutWOil_4 = QoilA_tau(i);
            Q_oilOutDispSep_4 = QoilC_tau(i);
        elseif i==5
            Q_oilOutWOil_5 = QoilA_tau(i);
            Q_oilOutDispSep_5 = QoilC_tau(i);
        end
 end
end


Q_oilOutDisp_tot=Q_oilOutDispSep_1+Q_oilOutDispSep_2+Q_oilOutDispSep_3+Q_oilOutDispSep_4+Q_oilOutDispSep_5;
Q_oilOutWoil_tot=Q_oilOutWOil_1+Q_oilOutWOil_2+Q_oilOutWOil_3+Q_oilOutWOil_4+Q_oilOutWOil_5;

% Estimating the waterflow out
Q_water_out = u(11)*rho_waterS/rho_water + u(2)*rho_oilS/rho_oil*wc + Q_oilOutWoil_tot; % Sm3/s -> m3/s

%the fractions are of the total flow going out of the separator and into
%the hydrocyclone
x_oilOutDispSep_1=Q_oilOutDispSep_1/Q_water_out; % Deviding on the total waterflow out, that is waterflow including dispersed oil
x_oilOutDispSep_2=Q_oilOutDispSep_2/Q_water_out;
x_oilOutDispSep_3=Q_oilOutDispSep_3/Q_water_out;
x_oilOutDispSep_4=Q_oilOutDispSep_4/Q_water_out;
x_oilOutDispSep_5=Q_oilOutDispSep_5/Q_water_out;

x_oilOutDispSep_tot=x_oilOutDispSep_1+x_oilOutDispSep_2+x_oilOutDispSep_3+x_oilOutDispSep_4+x_oilOutDispSep_5;

% This is the total oil fraction out with the water of the separator


%         Qoil_out_tot=Qoil_out0+Qoil_out+Q_oil_in;
%         Qoil_disp_out_tot=Qoil_disp_out0+Qoil_disp_out;

function WaterVol = WaterVol(R,Lweir,L,n,hweir,hw)
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 = -1;
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);


function f = hoBelowWeir(ho,R,Ldiff,n,Vo)
f = cyl(ho,R,Ldiff)+one_ellip(ho,R,n)-Vo;
f = (f/Vo);


function f = hoAboveWeir(ho,R,L,n,Vo,Vw,hweir,Lweir,VoWeir)
f = cyl(ho,R,L)+2*one_ellip(ho,R,n)-Vo-Vw;
f = (f/(Vo+Vw));


function f = hwBelowWeir(hw,R,Lweir,n,Vw)
f = cyl(hw,R,Lweir)+(pi/(2*n*n))*(R*hw*hw-(hw*hw*hw)/3)-Vw;
f = (f/Vw);


function f = hwAboveWeir(hw,R,L,n,Vw,VoWeir)
f = cyl(hw,R,L)+2*one_ellip(hw,R,n)-VoWeir-Vw;
f = (f/(Vw+VoWeir));

% 
function tau = residencetime(R_sep,Lweir,L_sep,n,hweir,hw0,q,delta_t)

% Knowing the volume of the water phase one can find the residence time 
% from the equation V= delta_t ' the sum from i=j(now) to (j-tau) (which
% is the time the volume element q went into the tank) over qi.

% q in the function is given as a vector. That is it's nessesary to have
% matlab store a number of values going in. Until enough flow has gone
% through to equal the amount of the watervolume, that is one retention 
% time, the program will not be able to run. 
% the 

% delta_t=10;

tau_m=0;
S=0; %the sum of the qs make up the volume the element has gone through
%When this volime equals the volume of the water phase it means that the 
%volume has been completely exchanged. That is that the little volume
%element going in is on it's way out. Giving the residence time. 

%Calculating the watervolume of the tank
% Vw=90; %Volume of the water phase

Vw=WaterVol(R_sep,Lweir,L_sep,n,hweir,hw0);

Vw_m=Vw/delta_t;

while S <= Vw_m
    S=S+q(end-tau_m);
%     disp(S);
    
    tau_m=tau_m+1;
    
end
tau_m = tau_m - 1;
tau = delta_t*tau_m; % i virkelig tid

return

