function sys=two_fluid(t,x,u,data,mode)
%   Model file for two-phase flow based on the two fluid model
%
%   Written by Espen Storkaas Jan. 15th 2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<5
    %for linearization
    if nargin==4
        u=[u;data];
        clear data
    end
    global data mode Y
end
persistent T
%Assigning inputs
x=x(:);

alpha_L=x(data.I_mL)/data.rho_L;
rho_G=x(data.I_mG)./(1-alpha_L);
W_L=x(data.I_iL);
W_G=x(data.I_iG);

z=u(1);
W_L_in=u(2);
W_G_in=u(3);
P_0=u(4);

%Computing derived variables
P = data.Z*rho_G*1e3*data.R*data.T/data.Mw_G; 
[W_L_choke,W_G_choke]=Sachdeva_choke(P(end),P_0,rho_G(end),alpha_L(end),z,data);
[u_L,u_G]=velocity(W_L,W_G,W_L_choke,W_G_choke,W_L_in,W_G_in,alpha_L,rho_G,z,data);
W_L=[W_L_in/data.A;W_L;W_L_choke/data.A];
W_G=[W_G_in/data.A;W_G;W_G_choke/data.A];
[M_L,M_G]=inertia(alpha_L,rho_G,u_L,u_G,data);
[g_L,g_G]=gravity(alpha_L,rho_G,data);
[F_wL,F_wG,F_i]=friction(alpha_L,rho_G,u_L,u_G,data);
% u_G
% W_G

%Conservation equations
if isequal(mode,'derivatives')
    sys(data.I_mL)=(W_L(data.I_mL)-W_L(data.I_mL+1))./data.dx_m;
    sys(data.I_mG)=(W_G(data.I_mL)-W_G(data.I_mL+1))./data.dx_m;
    sys(data.I_iL)=(M_L(data.I_i)-M_L(data.I_i+1)...
        +mean([alpha_L(data.I_i) alpha_L(data.I_i+1)],2).*(P(data.I_i)...
        -P(data.I_i+1)))./data.dx_i-g_L-F_wL+F_i;
    sys(data.I_iG)=(M_G(data.I_i)-M_G(data.I_i+1)...
        +mean([1-alpha_L(data.I_i) 1-alpha_L(data.I_i+1)],2).*(P(data.I_i)...
        -P(data.I_i+1)))./data.dx_i-g_G-F_wG-F_i;
    sys=[sys'];
    if isempty(T)
        T=t;
    elseif t-T>60
        T=round(t);
        disp(round(t/60))
    end
    if nargin<5
        rho_T=alpha_L(end)*data.rho_L+(1-alpha_L(end))*rho_G(end);
    W=W_G_choke+W_L_choke;
    Q=W/rho_T;
        %Y=[P/1e5;alpha_L;u_L;u_G;W_L_choke;W_G_choke;rho_T;W;Q];
      
    end
elseif isequal(mode,'measurements')
    %sys=[P;alpha_L;W_L*data.A;W_G*data.A;W_L_choke;W_G_choke];
    rho_T=alpha_L(end)*data.rho_L+(1-alpha_L(end))*rho_G(end);
    W=W_G_choke+W_L_choke;
    Q=W/rho_T;
    %sys=[P/1e5;alpha_L;u_L;u_G;W_L_choke;W_G_choke;rho_T;W;Q];
    sys=[P(1)/1e5;P(8)/1e5;P(13)/1e5;rho_T;W;Q];
end