function [sys,x0,str,ts]=flow_split(t,x,u,flag,rho_fluidS,rho_fluid,Diam1,Diam2,L_1,L_2,fric1,fric2,dz1,dz2)
%   Simulink S-function computing the flow in two merging pipes given the 
%   total flow and the inlet pressures.
%   See the file flow_split_basis.tex in addition to some handwritten notes for the background.

%   INPUTS:
%   U(1) Total flow in Sm3/s
%   U(2) Inlet pressure pipe 1 in Barg
%   U(3) Inlet pressure pipe 2 in Barg
%
%   STATES:
%   No states. (This is a static function.)
%
%   OUTPUTS:
%   Y(1) Flow rate in pipe 1 in Sm3/s
%   Y(2) Flow rate in pipe 2 in Sm3/s
%   Y(3) Pressure in merging point in Barg
%
%   PARAMETERS:
%   rho_fluidS  Standard fluid density  [kg/Sm3]
%   rho_fluid   Fluid density at pipe conditions [kg/m3]
%   Diam1       Diameter of pipe 1 [m]
%   Diam2       Diameter of pipe 2 [m]
%   L_1         Length of pipe 1 [m] 
%   L_2         Length of pipe 2 [m]
%   fric1       Friction factor pipe 1 [-]
%   fric2       Friction factor pipe 2 [-]
%   dz1         Height difference of pipe 1 [m] (out pos - in pos)        
%   dz2         Height difference of pipe 2 [m] (out pos - in pos)

%   See sfuntmpl.m for a general S-function template.
    
%   Copyright Oct 2001 ABB AS.
%
%   Rev 1.0 Oct 5th 2001 Olav Slupphaug.
%

switch flag

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

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1
    sys = [];

  %%%%%%%%%%%%%%%%%%%%%%%%
  % Update and Terminate %
  %%%%%%%%%%%%%%%%%%%%%%%%
  case {2,9}
    sys = []; % do nothing

  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3
    sys = mdlOutputs(u,rho_fluidS,rho_fluid,Diam1,Diam2,L_1,L_2,fric1,fric2,dz1,dz2); 

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

% end limintm

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes(xi)

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

sys = simsizes(sizes);
str = [];
x0  = [];
ts  = [0 0];   % sample time: [period, offset]

% end mdlInitializeSizes

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

sys=[];

% end mdlDerivatives

%
%=============================================================================
% mdlOutputs
% Return the output vector for the S-function
%=============================================================================
%
function sys = mdlOutputs(u,rho_fluidS,rho_fluid,Diam1,Diam2,L_1,L_2,fric1,fric2,dz1,dz2)
%   INPUTS:
%   U(1) Total flow in Sm3/s
%   U(2) Inlet pressure pipe 1 in Barg
%   U(3) Inlet pressure pipe 2 in Barg
%
%   OUTPUTS:
%   Y(1) Flow rate in pipe 1 in Sm3/s
%   Y(2) Flow rate in pipe 2 in Sm3/s
%   Y(3) Pressure in merging point in Barg
%
%   PARAMETERS:
%   rho_fluidS  Standard fluid density  [kg/Sm3]
%   rho_fluid   Fluid density at pipe conditions [kg/m3]
%   Diam1       Diameter of pipe 1 [m]
%   Diam2       Diameter of pipe 2 [m]
%   L_1         Length of pipe 1 [m] 
%   L_2         Length of pipe 2 [m]
%   fric1       Friction factor pipe 1 [-]
%   fric2       Friction factor pipe 2 [-]
%   dz1         Height difference of pipe 1 [m] (out pos - in pos)        
%   dz2         Height difference of pipe 2 [m] (out pos - in pos)

%   See the file flow_split_basis.tex in addition to some handwritten notes for the background.

gr      =   9.81;

g       =   u(1)*rho_fluidS;
a       =   u(2)*1e5;
d       =   u(3)*1e5;
b       =   rho_fluid*gr*dz1;
e       =   rho_fluid*gr*dz2;
c       =   rho_fluid*fric1*L_1/(Diam1*2)*(1/(rho_fluid*pi*(Diam1/2)^2))^2;
f       =   rho_fluid*fric2*L_2/(Diam2*2)*(1/(rho_fluid*pi*(Diam2/2)^2))^2;

% Start out by assuming positive flow direction in both pipes
root1    =   1/(2*(c-f))*(2*c*g+2*sqrt(e*c+c*a-c*b-c*d-e*f-f*a+f*b+d*f+f*c*g^2));
root2    =   1/(2*(c-f))*(2*c*g-2*sqrt(e*c+c*a-c*b-c*d-e*f-f*a+f*b+d*f+f*c*g^2));
z=0;x=0;
if isreal(root1),
    if min(root1,root2)<0,
       root=max(root1,root2); % Cannot choose the smaller one
    else
       root=min(root1,root2); % Choosing the larger one will probably yield x<0. This remains to be proved though.
    end
    z        =   root; 
    x        =   -root + g;
    y        =   -(-c*d+e*c+f*a-f*b-f*c*g^2+2*f*c*root*g)/(c-f);
	
    sys(1)  =   x/rho_fluidS;
    sys(2)  =   z/rho_fluidS;
    sys(3)  =   y/1e5;
end
if ~exist('sys') | z<0 | x<0, % meaning that the flow direction is not positive in both pipes. It is probably not 
                            % necessary to test on z<0 or x<0 because they will probably either be both positive 
                            % or both zero. This remains to be proved though.
    % Assume negative flow direction in pipe 1 and positive flow direction in pipe 2
    %$Z=\frac{1}{2\left( f+c\right) }\left( 2cg+2\sqrt{\left( -ef-fa+fb+df-fcg^{2}-ec-ca+cb+dc\right) }\right) $
    root1    =   1/(2*(c+f))*(2*c*g+2*sqrt(-e*c-c*a+c*b+c*d-e*f-f*a+f*b+d*f-f*c*g^2));
    root2    =   1/(2*(c+f))*(2*c*g-2*sqrt(-e*c-c*a+c*b+c*d-e*f-f*a+f*b+d*f-f*c*g^2));
    if isreal(root1),
        root=max(root1,root2);  % Cannot choose the smaller one, because it may cause a positive flow
                                % direction in pipe 1 
        z        =   root; 
        x        =   -root + g;
        %y=$\frac{dc-ec+fa-fb+fcg^{2}-2fc\rho g}{f+c}$
        y        =   (c*d-e*c+f*a-f*b+f*c*g^2-2*f*c*root*g)/(c+f);
        
        sys(1)  =   x/rho_fluidS;
        sys(2)  =   z/rho_fluidS;
        sys(3)  =   y/1e5;
    end   
    if ~exist('sys') | z>0 | x<0, % meaning that the flow direction is not negative in pipe 1 and positive in pipe 2.
                                % It is probably not necessary to test on z>0 or x<0 because they will probably 
                                % either satisfy z<=0 and x>=0 or both be zero. This remains to be proved though.
            % Assume positive flow direction in pipe 1 and negative flow direction in pipe 2
            % $Z=\frac{1}{2\left( f+c\right) }\left( 2cg+2\sqrt{\left( ef+fa-fb-df-fcg^{2}+ec+ca-cb-dc\right) }\right) $
            root1    =   1/(2*(c+f))*(2*c*g+2*sqrt(e*c+c*a-c*b-c*d+e*f+f*a-f*b-d*f-f*c*g^2));
            root2    =   1/(2*(c+f))*(2*c*g-2*sqrt(e*c+c*a-c*b-c*d+e*f+f*a-f*b-d*f-f*c*g^2));
            if isreal(root1),
                root=min(root1,root2);  % Cannot choose the larger one, because it may cause 
                                        % a positive flow direction in pipe 2
                z        =   root; 
                x        =   -root + g;
                %$y=\frac{dc-ec+fa-fb-fcg^{2}+2fc\rho g}{f+c}$
                y        =   (c*d-e*c+f*a-f*b-f*c*g^2+2*f*c*root*g)/(c+f);
                
                sys(1)  =   x/rho_fluidS;
                sys(2)  =   z/rho_fluidS;
                sys(3)  =   y/1e5;
            end  
    end                      
end
    
% end mdlOutputs

