function [sys,x0,str,ts] = sHydrocyclone(t,x,u,flag,rho_waterS,rho_water, ...
     rho_oil,rho_oilS,my,Dc,Dp,SC_underflow,n_liners)

% Inputs:
% u(1) - Pin, Pressure at inlet [barg]= bar over the 1 atm pressure you 
%        normally have on the earths surface.
% u(2) - Qoil_oS, Overflow of oil in standard conditions [Sm3/s]
% u(3) - Qoil_uS, Underflow of oil in standars conditions [Sm3/s]
% u(4) - xoil_in, Fraction of the inflow that is oil
% u(5) - x_oilInHyd_1, Fraction of flow to compartment 1
% u(6) - x_oilInHyd_2, Fraction of flow to compartment 2
% u(7) - x_oilInHyd_3, Fraction of flow to compartment 3
% u(8) - x_oilInHyd_4, Fraction of flow to compartment 4
% 
% Compartment 5, x_oilInHyd_5 will be given by the 4 other componenet 
% balances and the total mass balance Qoil_in = Qoil_u+Qoil_o
% 
% Outputs:
% sys(1) - Po, Overflow pressure [barg]
% sys(2) - Pu, Underflow pressure [barg]
% sys(3) - Qin, The inflow [Sm3/s]
% sys(4) -  N, hydrocyclone efficiency
% 
% Fractions of oil out with the overflow
% sys(5) - x_oilOutOverHyd_1, Fraction of flow overflow compartment 1
% sys(6) - x_oilOutOverHyd_2, Fraction of flow overflow compartment 2
% sys(7) - x_oilOutOverHyd_3, Fraction of flow overflow compartment 3
% sys(8) - x_oilOutOverHyd_4, Fraction of flow overflow compartment 4
% sys(9) - x_oilOutOverHyd_5, Fraction of flow overflow compartment 5
% 
% Fractions of oil out with the underflow
% sys(10) - x_oilOutUnderHyd_1, Fraction of flow overflow compartment 1
% sys(11) - x_oilOutUnderHyd_2, Fraction of flow overflow compartment 2
% sys(12) - x_oilOutUnderHyd_3, Fraction of flow overflow compartment 3
% sys(13) - x_oilOutUnderHyd_4, Fraction of flow overflow compartment 4
% sys(14) - x_oilOutUnderHyd_5, Fraction of flow overflow compartment 5
% 
% Parameters:
% rho_waterS          Water density at standard conditions [kg/Sm3]
% rho_water           Water density [kg/m3]
% rho_oilS            Oil density at standard conditions [kg/Sm3]
% rho_oil             Oil density [kg/m3]
% my                  Water viscosity
% Dc                  Hydrocyclone diameter at the broadest [m]
% Dp                  Particle size
% SC_underflow        Speed correction factor for the underflow. Correcting
%                     the diameter of the underflow being alot smaller than
%                     that of the cyclone [fraction 0-1]
% n_liners            number of liners in the cyclone
% 
%
% The following outlines the general structure of an S-function.
%
switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes;

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u);

  %%%%%%%%%%
  % Update %
  %%%%%%%%%%
  case 2,
    sys=mdlUpdate(t,x,u);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u,flag,rho_waterS,rho_water,rho_oil,rho_oilS,my,Dc, ... 
        Dp,SC_underflow,n_liners);

  %%%%%%%%%%%%%%%%%%%%%%%
  % GetTimeOfNextVarHit %
  %%%%%%%%%%%%%%%%%%%%%%%
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case 9,
    sys=mdlTerminate(t,x,u);

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

end

% end sfuntmpl

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts]=mdlInitializeSizes
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;

sizes.NumContStates  = 0;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 14; 
sizes.NumInputs      = 8; 
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

% initialize the initial conditions
%
x0  = [];

%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];

% end mdlInitializeSizes

%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)

sys = [];

% end mdlDerivatives

%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)

sys = [];

% end mdlUpdate

%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================


function sys=mdlOutputs(t,x,u,flag,rho_waterS,rho_water,rho_oil,rho_oilS,my, ...
     Dc,Dp,SC_underflow,n_liners)

%Initial Qin


Pin = u(1); %barg
QoS = u(2);
QuS = u(3);
% x_oil_in=u(4); %0-1
%QinS=u(2); %Sm3/s
%Qoil_inS=u(3); %Sm3/s
xoil_in = u(4);
x_oilInHyd_1 = u(5);
x_oilInHyd_2 = u(6);
x_oilInHyd_3 = u(7);
x_oilInHyd_4 = u(8);

x_oilInHyd_5 = xoil_in-(x_oilInHyd_1+x_oilInHyd_2+x_oilInHyd_3+x_oilInHyd_4);

Qo = QoS*rho_oilS/rho_oil/n_liners;
Qu = QuS*rho_oilS/rho_oil/n_liners;

% 1) Finding Qin from equation 1 knowing the flows going out and using a
% simple mass balance

Qin = Qo+Qu;

% 2a) Finding Qoil_in from equation 11 knowing the inflow as well as the 
% oil fraction of the inlet xoil_in from equation 

Qoil_in=xoil_in*Qin;

% 2b) Finding the Chamberflows
Qoil_in1=x_oilInHyd_1*Qin;
Qoil_in2=x_oilInHyd_2*Qin;
Qoil_in3=x_oilInHyd_3*Qin;
Qoil_in4=x_oilInHyd_4*Qin;
Qoil_in5=x_oilInHyd_5*Qin;

% 3) Finding deltaPo from equation 5
a=(rho_water/my)^2;
a2=a*Qin;
a3=a2/sqrt(rho_water);
b=Dc*sqrt(rho_water);
b2=b/my;
b3=0.2*(b2^1.85);
c=(a3/b3);
deltaPo=nthroot(c,0.425);

% 4) Finding Qwater_in from equation 8
Qwater_in=Qin-Qoil_in; % Transfer this equation to the separator eventually

% 5) Finding Po from equation 4
Po=Pin-(deltaPo/1e5); %bar

% 6) Finding Dpc from equation 3
d=abs(rho_water-rho_oil);
d2=d/(my^2);
e=Dc*sqrt(rho_water);
e2=(e/my)^0.85;
e3=0.177*e2;
Dpcm=sqrt(e3/(d2*(deltaPo^0.575)));

%%%

Qoil_inx=[Qoil_in1 Qoil_in2 Qoil_in3 Qoil_in4 Qoil_in5];
Dpm=Dp*1e-6; %m

Qoil_o1=0;
Qoil_o2=0;
Qoil_o3=0;
Qoil_o4=0;
Qoil_o5=0;

k=Dpm/Dpcm;
 
% 7) Finding the hydrocyclone efficiency from equation 9,
%    that is dependant on the cut sizes relationship to the size of the
%    particles

for i=1:length(Qoil_inx) 
     N(i)=(5/6)*k(i)-(1/3); % gir N=0 for 0.4 og N=1 ved k = 1.6
   if N(i)>1
       N(i) = 1;
   elseif N(i)<0
       N(i) = 0;
   else
       N(i)=N(i);
   end

% 8) Finding Qoil_o from equation 6
%    Qoil_o=N*Qoil_in;
   Qoil_o(i)=N(i)*Qoil_inx(i);

   
     % Keeping track of the compartment
     if i == 1
         Qoil_o1 = Qoil_o(i);

     elseif i == 2
         Qoil_o2 = Qoil_o(i);

     elseif i== 3
         Qoil_o3 = Qoil_o(i);

     elseif i == 4
         Qoil_o4 = Qoil_o(i);

     elseif i == 5
         Qoil_o5 = Qoil_o(i);
      
     end

end

   FS=Qo/Qin*100;
   if FS<1.53
   N_FS=59.5*FS/100;
    disp('FS too small')
   else
   N_FS=91/100;
   end
     
% 9) Finding Qoil_u from equation 7
%    Qoil_u=Qoil_in-Qoil_o;

Qoil_o1_FS = N_FS*Qoil_o1;
Qoil_u1_FS = Qoil_in1-Qoil_o1_FS;

Qoil_o2_FS = N_FS*Qoil_o2;
Qoil_u2_FS = Qoil_in2-Qoil_o2_FS;

Qoil_o3_FS = N_FS*Qoil_o3;
Qoil_u3_FS = Qoil_in3-Qoil_o3_FS;

Qoil_o4_FS = N_FS*Qoil_o4;
Qoil_u4_FS = Qoil_in4-Qoil_o4_FS;

Qoil_o5_FS = N_FS*Qoil_o5;
Qoil_u5_FS = Qoil_in5-Qoil_o5_FS;

Qoil_o_tot = Qoil_o1_FS + Qoil_o2_FS + Qoil_o3_FS + Qoil_o4_FS + Qoil_o5_FS;

Qoil_u_tot = Qoil_u1_FS + Qoil_u2_FS + Qoil_u3_FS + Qoil_u4_FS + Qoil_u5_FS;

if Qoil_o_tot >= Qo
    disp('Qoil_o > Qo')
    Qoil_o = Qo; % Qoil_o can't be larger than Qo
    Qoil_u = Qoil_u_tot+(Qoil_o_tot-Qoil_o);
    persentage_off = (Qoil_o-Qo)/Qo;
    
          Qoil_o1_FS = Qoil_o1_FS*(1-persentage_off);
          Qoil_o2_FS = Qoil_o2_FS*(1-persentage_off);
          Qoil_o3_FS = Qoil_o3_FS*(1-persentage_off);
          Qoil_o4_FS = Qoil_o4_FS*(1-persentage_off);
          Qoil_o5_FS = Qoil_o5_FS*(1-persentage_off);
          
          Qoil_u1_FS = Qoil_inx(i)-Qoil_o1_FS;
          Qoil_u2_FS = Qoil_inx(i)-Qoil_o2_FS;
          Qoil_u3_FS = Qoil_inx(i)-Qoil_o3_FS;
          Qoil_u4_FS = Qoil_inx(i)-Qoil_o4_FS;
          Qoil_u5_FS = Qoil_inx(i)-Qoil_o5_FS;
    
elseif Qoil_u_tot >= Qu
    disp('Qoil_u > Qu')
    Qoil_u = Qu; % Qoil_u cannot be larger than Qu
    Qoil_o = Qoil_o_tot+(Qoil_u_tot-Qoil_u);
    persentage_off = (Qoil_o-Qo)/Qo;
    
          Qoil_u1_FS = Qoil_u1_FS*(1-persentage_off);
          Qoil_u2_FS = Qoil_u2_FS*(1-persentage_off);
          Qoil_u3_FS = Qoil_u3_FS*(1-persentage_off);
          Qoil_u4_FS = Qoil_u4_FS*(1-persentage_off);
          Qoil_u5_FS = Qoil_u5_FS*(1-persentage_off);
          
          Qoil_o1_FS = Qoil_inx(i)-Qoil_u1_FS;
          Qoil_o2_FS = Qoil_inx(i)-Qoil_u2_FS;
          Qoil_o3_FS = Qoil_inx(i)-Qoil_u3_FS;
          Qoil_o4_FS = Qoil_inx(i)-Qoil_u4_FS;
          Qoil_o5_FS = Qoil_inx(i)-Qoil_u5_FS;
    
else
    Qoil_u = Qoil_u_tot;
    Qoil_o = Qoil_o_tot;
    % and the compartment flows remain as the originals
end

% 10) Finding the water flow understream, Qwater_u from equation 13
Qwater_u = Qu - Qoil_u;

% 11) Finding the waterfraction of the underflow, xwater_u  from equation 10
if Qu>0
xwater_u=Qwater_u/Qu; %eventually (1-xoil_u)

% 12) Finding the oilfraction of the underflow, xoil_u from equation 14.
xoil_u=Qoil_u/Qu;
else
   xwater_u=0; %eventually (1-xoil_u)

% 12) Finding the oilfraction of the underflow, xoil_u from equation 14.
    xoil_u=0; 
end
% 13)Finding the underflow density rho_u from 16
rho_u=(xoil_u*rho_oil)+(xwater_u*rho_water);

% 14) Finding the oilfraction of the overflow from equation 17
xoil_o=Qoil_o/max(Qo,1e-8);

% 15) Finding the water fraction at the inlet, xwater_in from 12
if Qin>0
xwater_in=Qwater_in/Qin;
else
    xwater_in=0;
end

% 16) Finding the inlet density rho_in from 15 
rho_in=(xoil_in*rho_oil)+(xwater_in*rho_water);


% The height of the inflow is then the hight of half a hydrocyclone radius.
% as it is placed horizontaly. 

% The height of the underflow is zero as the reference frame has been
% placed here. 

hin=0.5*Dc;%10*Dc-((0.26/2)*Dc);

g = 9.81; %The gravity constant [m/s2]

% 17) Finding Pu from equation 2 Bernoulli corrected 
Pu=rho_u*((((4*Qin/(pi*(Dc^2)))^2)/2)-(((4*Qu/(pi*((Dc/SC_underflow)^2)))^2)/2)+(g*hin)+(Pin*1e5/rho_in)); % the pressure in Pa for the right denomination
% ((Qu*300)^2))/2)
% Pu=(((Qin^2)/2)+g*hin+(Pin*1e5/rho_in)-(((xoil_u*Qu)^2)/2))*rho_u; %Pa
deltaPuPa = (Pin*1e5)-Pu;
deltaPu = deltaPuPa/1e5;  %[barg]


if Qoil_in>0
Ntot= Qoil_o/Qoil_in;
else
    Ntot=0;
end

if Qo>0
    x_oilOutOverHyd_1 = Qoil_o1_FS/Qo;
    x_oilOutOverHyd_2 = Qoil_o2_FS/Qo;
    x_oilOutOverHyd_3 = Qoil_o3_FS/Qo;
    x_oilOutOverHyd_4 = Qoil_o4_FS/Qo;
    x_oilOutOverHyd_5 = Qoil_o5_FS/Qo;
else
    x_oilOutOverHyd_1 = 0;
    x_oilOutOverHyd_2 = 0;
    x_oilOutOverHyd_3 = 0;
    x_oilOutOverHyd_4 = 0;
    x_oilOutOverHyd_5 = 0;
end
if Qu>0
x_oilOutUnderHyd_1 = Qoil_u1_FS/Qu;
x_oilOutUnderHyd_2 = Qoil_u2_FS/Qu;
x_oilOutUnderHyd_3 = Qoil_u3_FS/Qu;
x_oilOutUnderHyd_4 = Qoil_u4_FS/Qu;
x_oilOutUnderHyd_5 = Qoil_u5_FS/Qu;
else
x_oilOutUnderHyd_1 = 0;
x_oilOutUnderHyd_2 = 0;
x_oilOutUnderHyd_3 = 0;
x_oilOutUnderHyd_4 = 0;
x_oilOutUnderHyd_5 = 0;
end

% Outputs
sys(1)=Po;
sys(2)=Pu/1e5;
sys(3)=Qin*rho_oil*n_liners/rho_oilS; %[Sm3/s];
%sys(4)=Qwater_out*rho_oil/rho_oilS;
% sys(4)=xoil_u;
% Oil flows out with the overflow
sys(4)=Ntot;
sys(5)=x_oilOutOverHyd_1;
sys(6)=x_oilOutOverHyd_2;
sys(7)=x_oilOutOverHyd_3;
sys(8)=x_oilOutOverHyd_4;
sys(9)=x_oilOutOverHyd_5;

% Oil flow out with the underflow
sys(10)=x_oilOutUnderHyd_1; 
sys(11)=x_oilOutUnderHyd_2; 
sys(12)=x_oilOutUnderHyd_3; 
sys(13)=x_oilOutUnderHyd_4; 
sys(14)=x_oilOutUnderHyd_5; 


% end mdlOutputs

% 
% =============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block.  Note that the result is
% absolute time.  Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
% =============================================================================
% %
function sys=mdlGetTimeOfNextVarHit(t,x,u)

sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

% end mdlGetTimeOfNextVarHit

%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

    
