function [sys]=slug_3D(t,x,u,output,data)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Model file for the simplified slug model
%   For SIMULINK use, called by s_slug.m
%   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}]
%   
%
%   Written by :    Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
%   Date:           08.10.02
%   
%   For questions, comments or documentation, contact author
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 




%Assingning states and inputs
mL=x(1);
mG1=x(2);
mG2=x(3);
%Inputs
z=u(1);
mG_in=u(2);
P_feed=u(3);
P0=u(4)*1e5;
%initializing implicit rutine


%Model

P1=mG1*data.R*data.T/(data.M_G*data.V_G1);
mL_in=data.K4*sqrt((P_feed>P1)*(P_feed-P1));
rho_G1=mG1/data.V_G1;
V_L=mL/data.rho_L;
%Second order equation to compute rho_mix
beta=data.g*(data.H2+data.H3);
a=-beta*data.V_T/data.rho_L-data.g*data.V_T^2/(data.A1*data.rho_L);
b=data.V_T*(beta+V_L*data.g/data.A1+data.g*data.V_T/data.A1+2*data.g*mG2/(data.rho_L*data.A1)+P1/data.rho_L)+beta/data.rho_L*mG2;
c=-data.rho_L*data.g*V_L*data.V_T/data.A1-data.g*(V_L+data.V_T)*mG2/data.A1-...
    data.g*mG2^2/(data.A1*data.rho_L)-P1*data.V_T-P1*mG2/data.rho_L+mG2*data.R*data.T/data.M_G;
rho_mix=real(min(roots([a b c])));

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);
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>0)*(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);
y=P0/P2;x1=1-alpha_Lm;
mdot_mix=(y<1)*data.K1*z*sqrt(P2*rho_T^2*...
    ((1-x1)*(1-y)/data.rho_L+x1*data.k*(1-y^((data.k-1)/data.k))/((data.k-1)*rho_G2)));
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;

if isequal(output,'derivatives')
    sys=[dmLdt;dmG1dt;dmG2dt];%t
    if ~isreal(sys)
        disp('Complex')
        sys=0*sys;
    end
elseif isequal(output,'measurements')
    sys=[P1/1e5 (P2-P0)/1e5 rho_T mdot_mix mdot_mix/rho_T];
        if ~isreal(sys)
        disp('Complex')
        sys=0*sys;
    end
end 
