function [x_dot,y] = v1_gaslift_model0(x,u,par)
%% **** States **********
% Check for unpyshical state input
for i=1:1:length(x)
    if x(i) <0
        x(i)=0;
    end
end
% Check for unreal state input
if( ~isreal(x))
    x=0*x;
end

x1=x(1);        		%Gas mass in annulus
x2=x(2);        		%Gas mass in riser
x3=x(3);        		%Liquid mass in riser (whole system, tubing+BH)	


%% ***** Inputs ******
u1=u(1);    			%Valve opening top choke
u2=u(2);    			%Valve opening of to injection choke
%% ****** Disturbances ******
GOR = u(3);        % Mass GOR from reservoir (d1) [-]
P_r = u(4)*1e5;    % Rservoir pressure (d2) [bar]
P_gs = u(5)*1e5;   % Gass source pressure (d3) [bar]
P0 = u(6)*1e5;     % Separator Pressure (d4) [bar] 

%% Density of gas inside riser
%rho_G_r=x2/(par.V_r + par.L_bh*par.S_bh - x3/par.rho_L);
if (par.V_r - (x3-par.rho_L*par.L_bh*par.S_bh)/par.rho_L) >0
    rho_G_r=x2/(par.V_r - (x3-par.rho_L*par.L_bh*par.S_bh)/par.rho_L);
else
    rho_G_r=0;
end

%% ************* Average mixture density ***************
%rho_mix=(x2+x3-par.rho_L*par.L_bh*par.S_bh)/(par.V_r);
rho_mix=(x2+x3-par.rho_L*par.V_bh)/(par.L_r*pi*par.r_r^2); rho_mix=max(rho_mix,0);
alpha_L_av = x3/(par.rho_L*(par.V_bh+par.V_r)); alpha_L_av=max(alpha_L_av,0);
%% Pressure riser top   
P_r_t=rho_G_r*par.R*par.T_r/(par.M_G_r_t); P_r_t=max(P_r_t,0);
%P_r_t = max(P0,P_r_t);

alpha_G_m_in = GOR/(GOR+1);
rho_G_r_in = (P_r - par.rho_L*par.g*par.L_bh)*par.M_G_r_t/(par.R*par.T_r);
%% ********* Pressure at bottom of annulus *********************
P_a_t = (par.R*par.T_a*x1/(par.M_G_a*par.V_a));
P_a_b= P_a_t + (x1*par.g*par.L_a/par.V_a) ; 

%% ********* Density of gas in the gas injection tank ************
rho_G_in= P_gs*par.M_G_a/(par.R*par.T_a);

%% ********** Injected gas into annulus ******************
%w_G_a_in = par.K_s*u2*sqrt(rho_G_in*max(P_gs-P_a_t,0));
w_G_a_in =u2; %Replace in this model beacuse we use constat gas into annulus. Normally use the one over
%% ********** Density of gas in bottom of annulus **********
rho_G_a_b= P_a_b*par.M_G_a/(par.R*par.T_a); 

%% ************  Calculating Friction of Riser **********
w_r_in_av = 14.6;
Usl_av=(1-alpha_G_m_in)*w_r_in_av/(par.rho_L*pi*(par.r_r^2)); Usl_av=max(Usl_av,0);
Usg_av= (alpha_G_m_in*w_r_in_av+w_G_a_in)/(rho_G_r*pi*(par.r_r^2)); Usg_av=max(Usg_av,0);
U_M_av=Usl_av+Usg_av; 
Re=(2*rho_mix*U_M_av*par.r_r)/par.my; Re=max(Re,0);
%lambda=0.0056+(0.5*(Re^(-0.32)));
%% Darcy Friction Factor
temp = -1.8*log10((par.ew/par.D_w/3.7)^1.11+6.9/Re);
lambda=(1/temp)^2;
%F_riser=(alpha_L_r_avg*lambda*rho_mix*(U_M_av^2)*(par.L_r+par.L_h))/(4*par.r_r);
F_riser=real(alpha_L_av*(lambda*rho_mix*(U_M_av^2)*(par.L_r+par.L_h))/(4*par.r_r));
F_riser=max(F_riser,0);

%% *********** Pressure at bottom of riser ******************
P_r_b=P_r_t + (rho_mix*par.g*par.L_r) +F_riser;  
if P_r_b<0 %Check for positive pressure in bottom of riser
   disp('Error: Pressure in bottom of riser is negative');
   P_r_b=0;
end


%% *********  Injection gas flow rate from annulus to riser	*********
if (P_a_b - P_r_b)>0 && rho_G_a_b>0 
    w_G_a_out = par.K_a*sqrt(rho_G_a_b*(P_a_b - P_r_b)); 
    %ORF = abs(1/(1^2*par.CD^2) -1);
    %w_G_a_out = par.A_v*sqrt(2*rho_G_a_b*max(0,P_a_b - P_r_b)/ORF);
else
    %disp('Error: Injection gas from annulus has negative pressure difference under root sign');
    w_G_a_out = 0; %Gas mass flow out of annulus
end

%% Liquid velocity at bottom hole
Ul_b= w_r_in_av/(par.rho_L*pi*(par.r_r^2)); Ul_b=max(Ul_b,0);
%% Reynolds number at bottom-hole:
Re_b=(2*par.rho_L*Ul_b*par.r_r)/par.my; Re_b=max(Re_b,0);
%% Friction factor at bottom-hole:
%lambda_b=0.0056+(0.5*(Re_b^(-0.32)));
%% Darcy Friction Factor
temp = -1.8*log10((par.ew/par.D_b/3.7)^1.11+6.9/Re_b);
lambda_b=(1/temp)^2;
%% Pressure loss due to friction from injection point to bottom-hole:
F_b=real((lambda_b*par.rho_L*(Ul_b^2)*(par.L_bh))/(4*par.r_r));
F_b=max(0,F_b);
%% Bottom-hole pressure
P_bh=P_r_b+par.rho_L*par.g*par.L_bh +F_b; P_bh=max(P_bh,0);

%% Liquid Inflow rate
if (P_r_b >= P_r)
    w_r_in=0;
else
    w_r_in=par.PI*(P_r-P_bh);
end
w_L_r_in = (1-alpha_G_m_in)*w_r_in; w_L_r_in=max(w_L_r_in,0);
w_G_r_in = alpha_G_m_in*w_r_in; w_G_r_in=max(w_G_r_in,0);
%% **********  Gas density in bottom of riser **************
rho_G_r_b=P_r_b*par.M_G_r_t/(par.R*par.T_r); 
rho_G_r_b=max(rho_G_r_b,0);

%% Alpha liquid in 
%alpha_G = x2/(x2+x3);
    
if ( ((w_G_a_out+w_G_r_in) >0)&&(w_L_r_in >0) )
    alpha_L_in= (w_L_r_in*rho_G_r_b)...  						
        /(w_L_r_in*rho_G_r_b+(w_G_a_out+w_G_r_in)*par.rho_L) ;
elseif (w_L_r_in >0)
    alpha_L_in=1;												
else
    alpha_L_in=0;
end
alpha_L_in=max(alpha_L_in,0);
alpha_L_in=min(alpha_L_in,1);
%% ******** Liquid volume fraction top of riser ****************       
if ( (2*alpha_L_av - alpha_L_in )> 0 )
    alpha_L_t=2*alpha_L_av - alpha_L_in; 	
else
    %disp('Error: Liquid volume fraction on top of riser is faulty calculated');
    alpha_L_t=0;
end

alpha_L_t=min(1,alpha_L_t); %Because the approximation over sometimes gives higher value than 1
alpha_L_t=max(0,alpha_L_t);
%% ***********  Density mixture top of riser *************************
rho_M_t= alpha_L_t*par.rho_L +(1-alpha_L_t)*rho_G_r; 	
rho_M_t=max(0,rho_M_t);
%% ********* Liquid mass fraction top of riser ********************
if ((alpha_L_t*par.rho_L +(1-alpha_L_t)*rho_G_r) > 0)
    alpha_L_m_t=(alpha_L_t*par.rho_L)...   			
        /(alpha_L_t*par.rho_L +(1-alpha_L_t)*rho_G_r);
else
    %disp('Error: Liquid mass fraction on top of riser is faulty calculated');
    alpha_L_m_t=0;
end
alpha_L_m_t=min(1,alpha_L_m_t);
alpha_L_m_t=max(0,alpha_L_m_t);
alpha_G_m_t=1-alpha_L_m_t;

%% **********  Mixture flow out of riser ****************************
if (P_r_t-P0)>0
    if rho_M_t<0
        disp('Error: Mixture Density on top of riser is lesser than zero')
    end
    if(rho_M_t==0)
        w_mix_r_out=0;
        disp('Error: Mixture Density on top of riser is zero')
    else
        w_mix_r_out=par.K_r*u1*sqrt(rho_M_t*(P_r_t-P0)); 	
        %ORF = abs(1/(u1^2*par.CD^2) -1);
        %w_mix_r_out = par.A_c*sqrt(2*rho_M_t*max(0,P_r_t-P0)/ORF);
    end   
else
    w_mix_r_out=0; 					
end
w_mix_r_out=max(0,w_mix_r_out);
%% Mass liquid flow rate out of riser
w_L_r_out=(1-alpha_G_m_t)*w_mix_r_out;
w_L_r_out=max(w_L_r_out,0);
w_L_r_out=min(w_L_r_out,w_mix_r_out);
%% Mass gas flow rate out of riser
w_G_r_out=alpha_G_m_t*w_mix_r_out;
w_G_r_out=max(0,w_G_r_out);

%% Derivatives
x1_dot = w_G_a_in - w_G_a_out;              							
x2_dot = (w_G_a_out+w_G_r_in) - w_G_r_out;  							
x3_dot = w_L_r_in - w_L_r_out;                                        

%% Derivatives vector
x_dot=[x1_dot;x2_dot;x3_dot]; 											

%% Measurements vector
%y=[P_r_b; P_r_t; P_a_b; w_mix_r_out; Q_out;w_G_a_out;w_L_r_in;alpha_G;alpha_G_m_t];
y=[w_G_a_in;P_a_t;P_r_t;P_a_b;P_bh;w_mix_r_out;rho_M_t;alpha_L_t];

%% ****** Checking of answers **************%
%Check for real states
if( ~isreal(x_dot))
    disp('error in model, x_dot not real, x_dot is: ');
    x_dot 
    x_dot=0*x_dot;
end
%Check for negative states (unphysical)
if(x<0)
    x_dot=0*x_dot;
    disp('error in model, x lesser than 0')
end
% Check for real measurments
if( ~isreal(y))
    disp('error in model, y not real')
    y=0*y;
end
%Check for Inf of NaN
for i=1:length(x_dot)
    if ((x_dot(i) == Inf) || isnan(x_dot(i)))
        x_dot(i)=0;
    end
end

    
end