function [x0,y,data]=initialize(z,mG_in,mL_in,P1,h1,P2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Initialization file for the simplified slug model
%   Performs a initialization to the stationary point for the given inputs
% 
%   Inputs: u=  z       Valve position
%               mG_in   Gas feed in kg/s
%               mL_in   Liquid feed in kg/s
%               P1      Upstream pressure
%               h1      Liquid level at bend
%               P2      Downstream pressure
%               
%               
%   Output: y=  x0      State vector
%               [rho_mix;v_G1;alpha_LT;V_LR;V_G2;mdot_mix;rho_T;v_G1] -  Various measurements 
%               data    struct containing case data
%
%   Written by :    Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
%   Date:           08.10.02
%   
%   For questions, comments or documentation, contact author
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options=optimset('fsolve');options=optimset('display','off','TolFun',1e-8,'TolX',1e-8);
%Generating Case data
gen_data
%equations
mdot_G1=mG_in;                                                         
mdot_Lout=mL_in;                                                        
mdot_Gout=mdot_G1;                                                      
alpha_Lm=1/(mG_in/mL_in+1);                                             
mdot_mix=mdot_Gout/(1-alpha_Lm);                                        
rho_mix=(P1-P2+data.rho_L*data.g*h1)/(data.g*(data.H2+data.H3));                                 %Eq. 11
V_LR=fsolve(@imp_fun,0.5*data.V_T,options,P2,rho_mix,data);     
V_G2=data.V_T-V_LR;                                                          
mG2=P2*V_G2*data.M_G/(data.R*data.T);                                                  
rho_G2=mG2/V_G2;                                                        
alpha_LT=alpha_Lm*rho_G2/((1-alpha_Lm)*data.rho_L+alpha_Lm*rho_G2);          
rho_T=alpha_LT*data.rho_L+(1-alpha_LT)*rho_G2;                               
alpha_L=V_LR/data.V_T;                                                       
mG1=P1*data.V_G1*data.M_G/(data.R*data.T);                                                  
rho_G1=mG1/data.V_G1;                                                        


%=========================================================================\
% 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));
    A1=data.r^2*(pi-phi-cos(pi-phi)*sin(pi-phi));  
    A=A1;
elseif(h1>=data.H1)
    A1=0;
    A=A1;
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=mdot_G1/(rho_G1*A);                                                 
V_L=V_LR+h1*data.A1;                                                         
mL=V_L*data.rho_L;                                                           
data.K3=1.1931%espen=1.1931; %fsolve(@imp_fun1,(data.rho_L-rho_G1)/(2*rho_G1*v_G1^2),options,... %paper=data.K3=0.11;
    %v_G1,rho_G1,V_LR,alpha_L,alpha_LT,data);      
data.K2=0.8619;%espen=0.8619; %v_G1/((data.H1-h1)/data.H1*sqrt((P1-P2-data.rho_L*data.g*data.H2*alpha_L)/rho_G1));    %paper=data.K2=6.84;       
data.K1=0.0050;%espen =0.0050; %mdot_mix/(z*sqrt(rho_T*(P2-data.P0))); %paper=data.K1=0.0054;
data.rho_mix_i=rho_mix;

%=========================================================================\
% heisann hoppsann fallereaal
%=========================================================================/


x0=[mL;mG1;mG2];
y=[rho_mix;v_G1;alpha_LT;V_LR;V_G2;mdot_mix;rho_T;v_G1];
if h1>data.H1
    disp('h1 too large, try again')
    h1
elseif h1<0
    disp('h1 negative, try again')
    h1
end

if (abs(mdot_Gout-mG_in)>0.001 | abs(mdot_Lout-mL_in)>0.001)
    disp('Something is wrong here!')
    mdot_Gout
    mdot_Lout

end

function diff=imp_fun(V_LR,P2,rho_mix,data)
V_G2=data.V_T-V_LR;                                                              
mG2=P2*V_G2*data.M_G/(data.R*data.T);                                                      
diff=(mG2+V_LR*data.rho_L)/data.V_T-rho_mix;                                         

function diff=imp_fun1(K3,v_G1,rho_G1,V_LR,alpha_L,alpha_LT,data);
w=(K3*rho_G1*v_G1^2/(data.rho_L-rho_G1))^data.n;
diff=((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)))-alpha_LT;