function [sys]=slug_3D_lin(t,x,u)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Model file for the simplified slug model
%   For linearization of model, ODE structure
% 
%   States: x=[mL;mG1;mG2]
%   Inputs: u=[z{-};mG_in{kg/s};mL_in{kg/s}]
%   Output: y=[P1{Bar};P2{Bar};rho_T{kg/m3};Q{m3/s};W{kg/s}]
%   Dependencies:   -Case data in global struct 'data' 
%                    (Can be generated with initialize.m)
%   
%   Linearization can be performed with linearization.m
%
%   Written by :    Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
%   Date:           08.10.02
%   
%   For questions, comments or documentation, contact author
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global Y data 
persistent rho_mix options

%Assingning states and inputs
mL=x(1);
mG1=x(2);
mG2=x(3);
%Inputs
z=u(1);
mG_in=u(2);
mL_in=u(3);
%initializing implicit rutine
if isempty(rho_mix)
    rho_mix=data.rho_mix_i;
end
if isempty(options)
    options=optimset('fsolve');options=optimset('display','off','TolFun',1e-10,'TolX',1e-10);
end

%Model

P1=mG1*data.R*data.T/(data.M_G*data.V_G1);
rho_G1=mG1/data.V_G1;
V_L=mL/data.rho_L;
rho_mix=fsolve(@imp_func,rho_mix,options,V_L,mG2,P1,data);
V_LR=(rho_mix*data.V_T-mG2)/data.rho_L;
h1=(V_L-V_LR)/data.A1;
V_G2=data.V_T-V_LR;
P2=mG2*data.R*data.T/(V_G2*data.M_G);
alpha_L=V_LR/data.V_T;
phi=((data.H1-h1)*cos(data.theta)<data.r)*(pi-acos(1-(data.H1-h1)*cos(data.theta)/data.r))+...
    ((data.H1-h1)*cos(data.theta)>data.r)*(acos((data.H1-h1)*cos(data.theta)/data.r-1));
A=data.r^2*(pi-phi-cos(pi-phi)*sin(pi-phi));
v_G1=(h1<data.H1)*data.K2*(data.H1-h1)/data.H1*sqrt((P1-P2-data.rho_L*data.g*data.H2*alpha_L)/rho_G1);
mdot_G1=v_G1*rho_G1*A;
rho_G2=mG2/V_G2;
w=(data.K3*rho_G1*v_G1^2/(data.rho_L-rho_G1))^data.n;
alpha_LT=(V_LR>data.A2*data.H2)*(V_LR-data.A2*data.H2)/(data.A3*data.H3)+w/(1+w)...
    *(alpha_L-(V_LR>data.A2*data.H2)*(V_LR-data.A2*data.H2)/(data.A3*data.H3));
rho_T=alpha_LT*data.rho_L+(1-alpha_LT)*rho_G2;
alpha_Lm=alpha_LT*data.rho_L/(alpha_LT*data.rho_L+(1-alpha_LT)*rho_G2);
mdot_mix=data.K1*z*sqrt(rho_T*(P2-data.P0));
mdot_Lout=alpha_Lm*mdot_mix;
mdot_Gout=(1-alpha_Lm)*mdot_mix;
%Conservation equations
dmG1dt=mG_in-mdot_G1;
dmLdt=mL_in-mdot_Lout;
dmG2dt=mdot_G1-mdot_Gout;


sys=[dmLdt;dmG1dt;dmG2dt];
Y=[P1/1e5 P2/1e5 rho_T mdot_mix mdot_mix/rho_T];

function [sys]=imp_func(rho_mix,V_L,mG2,P1,data)
V_LR=(rho_mix*data.V_T-mG2)/data.rho_L;
h1=max((V_L-V_LR)/data.A1,0);
V_G2=data.V_T-V_LR;
P2=mG2*data.R*data.T/(V_G2*data.M_G);
sys=rho_mix*data.g*(data.H2+data.H3)-data.rho_L*data.g*h1-P1+P2;



