function [PipelineData]=ComputeConstants(InitCond,PipelineData)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 {Bar}
%               h1      Liquid level at bend
%               P2      Downstream pressure {Bar}
%               P0      Recieving (separator) pressure {Bar}
%               
%               
%   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
%   Dependencies:   gen_data.m
%
%   Written by :    Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
%   Date:           08.10.02
%   
%   For questions, comments or documentation, contact author
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options=optimset('fzero');options=optimset('display','off','TolFun',1e-8,'TolX',1e-8);
% Generating Case data
% converting from bar to Pa
P0=InitCond.P3*1e5;P1=InitCond.P1*1e5;P2=InitCond.P2*1e5;
%equations
mdot_G1=InitCond.mG_in;                                                         
mdot_Lout=InitCond.mL_in;                                                        
mdot_Gout=mdot_G1;                                                      
alpha_Lm=1/(InitCond.mG_in/InitCond.mL_in+1);                                             
mdot_mix=mdot_Gout/(1-alpha_Lm);                                        
rho_mix=(P1-P2+PipelineData.rho_L*PipelineData.g*InitCond.h1)/(PipelineData.g*(PipelineData.H2+PipelineData.H3));                                 %Eq. 11
V_LR=fzero(@imp_fun,0.5*PipelineData.V_T,options,P2,rho_mix,PipelineData);     
V_G2=PipelineData.V_T-V_LR;                                                          
mG2=P2*V_G2*PipelineData.M_G/(PipelineData.R*PipelineData.T);                                                  
rho_G2=mG2/V_G2;                                                        
alpha_LT=alpha_Lm*rho_G2/((1-alpha_Lm)*PipelineData.rho_L+alpha_Lm*rho_G2);          
rho_T=alpha_LT*PipelineData.rho_L+(1-alpha_LT)*rho_G2;                               
alpha_L=V_LR/PipelineData.V_T;                                                       
mG1=P1*PipelineData.V_G1*PipelineData.M_G/(PipelineData.R*PipelineData.T);                                                  
rho_G1=mG1/PipelineData.V_G1;                                                        
phi=((PipelineData.H1-InitCond.h1)*cos(PipelineData.theta)<PipelineData.r)*(pi-acos(1-(PipelineData.H1-InitCond.h1)*cos(PipelineData.theta)/PipelineData.r))+...        
    ((PipelineData.H1-InitCond.h1)*cos(PipelineData.theta)>PipelineData.r)*(acos((PipelineData.H1-InitCond.h1)*cos(PipelineData.theta)/PipelineData.r-1));
A=PipelineData.r^2*(pi-phi-cos(pi-phi)*sin(pi-phi));                                 
v_G1=mdot_G1/(rho_G1*A);                                                 
V_L=V_LR+InitCond.h1*PipelineData.A1;                                                         
mL=V_L*PipelineData.rho_L;                                                           
PipelineData.K3=fzero(@imp_fun1,(PipelineData.rho_L-rho_G1)/(2*rho_G1*v_G1^2),options,...
    v_G1,rho_G1,V_LR,alpha_L,alpha_LT,PipelineData);      
PipelineData.K2=v_G1/((PipelineData.H1-InitCond.h1)/PipelineData.H1*sqrt((P1-P2-PipelineData.rho_L*PipelineData.g*PipelineData.H2*alpha_L)/rho_G1));           
PipelineData.K1=mdot_mix/(InitCond.z*sqrt(rho_T*(P2-P0)));
if PipelineData.PdepFeed
    PipelineData.K4=(InitCond.mL_in+InitCond.mG_in)/(InitCond.P_f-InitCond.P1);
end
if InitCond.h1>PipelineData.H1
    disp('InitCond.h1 too large, try again')
    InitCond.h1
elseif InitCond.h1<0
    disp('InitCond.h1 negative, try again')
    InitCond.h1
end

if (abs(mdot_Gout-InitCond.mG_in)>0.001 || abs(mdot_Lout-InitCond.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;




