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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
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;
%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=min(roots([a b c]));

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;


%=========================================================================\
% Opprinnelig måte å beregne arealet på, men endret til h1<data.H1 :
%=========================================================================/
if(h1<data.H1)
    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));  
elseif(h1>=data.H1)
    A=0;
end

%=========================================================================\
% Forenklet måte å beregne arealet på
%=========================================================================/
% A=(data.H1>h1)*(data.H1-h1)*2*data.r;

% %=========================================================================\
% % Forsøk på å forbedre måten å beregne arealet på : 
% %=========================================================================/
% D_bend=data.H1-data.r*tan(data.theta);
% R_bend=D_bend/2;
% if(h1 >= data.H1)
%     A3=0;%gassen blir blokkert, areal lik 0
% elseif(h1<data.H1 & h1>R_bend) %her ligger forenklingen om at  røret ligger flatt. (radius bør være R_bend)
%     c=2*sqrt(h1*(2*R_bend-h1));
%     if(~isreal(c))
%         %disp('c er kompleks, h1>R_Bend')
%     end
%     theta_bend=2*asin(c/(2*R_bend));
%      if(~isreal(theta_bend))
%         %disp('theta blir kompleks for h1<R_bend')
%     end
%     
%     A_trekant=c*(h1-R_bend)/2;
%     A3=0.5*R_bend^2*theta_bend - A_trekant;
% elseif(h1<data.H1 & h1==R_bend)
%     A3=(pi*R_bend^2)/2;
%     
% elseif(h1<data.H1 & h1<R_bend)
%     c=2*sqrt(h1*(2*R_bend-h1));
%     if(~isreal(c))
%         %disp('c er kompleks, h1<R_Bend')
%     end
%     theta_bend=2*asin(c/(2*R_bend));
%     
%     if(~isreal(theta_bend))
%         %disp('theta blir kompleks for h1<R_bend')
%     end
%     
%     A_trekant=c*(R_bend-h1)/2;
%     A3=0.5*R_bend^2*theta_bend-A_trekant;
% end
% 
% if(~isreal(A3)) %BRUTE FORCE!!!!!!!! HAR IKKE FUNNET HVORFOR MAN FÅR 0+0.000000001i innimellom
%     A3=0;
% end
% A=A3;

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;

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