%Track_stationary.m
clear all
global data IX IU IFIXU uorg
IX=1:4;
IU=5:7;
IFIXU=2:3;
%Starting Initaialization
mG_in=0.362;
mL_in=8.64;
z=0.13;
P1_stasj=71.7e5;
h1_stasj=0.1;
P2_stasj=53.5e5;
[x0,y_stasj,data]=initialize(z,mG_in,mL_in,P1_stasj,h1_stasj,P2_stasj);
%End Initialization
%2 starting points must be generated
%This is done by running slug_3D_script to stationary soultion and saving
%the last value as initial value
%Commands: (After simulation)
%   x0=x(:,end)
%   save xinit... x0

%Starting point 1
z=0.07;
load xinit007
u=[z;mG_in;mL_in];
z1=[x0';u];
%starting point 2
load xinit0075
z=0.075
u=[z;mG_in;mL_in];
uorg=u;
z2=[x0';u];
%
[out,uu,step1,tt]=secant_cont_circle('slug_wrapper',z2,z1,0.01);
disp('Computing measurements')
xs=out(:,IX)';us=out(:,IU)';
for i=1:length(out)
    y(i,:)=slug_3DAE(0,xs(:,i),us(:,i),'measurements',data);
end
load slugdata
load OLGAdata

figure(1)
plot(us(1,:),y(:,3),'r--',z_vec,maxP1,'r-',z_vec,minP1,'r-',zOLGA,P1maxOLGA*1e5,'b-',zOLGA,P1minOLGA*1e5,'b-',zOLGA,P1stasjOLGA*1e5,'b--')
figure(2)
plot(us(1,:),y(:,4),'r--',z_vec,maxP2,'r-',z_vec,minP2,'r-',zOLGA,P2maxOLGA*1e5,'b-',zOLGA,P2minOLGA*1e5,'b-',zOLGA,P2stasjOLGA*1e5,'b--')
P1=y(:,3);
P2=y(:,4);
z=us(1,:);
rho_mix=xs(4,:);
save Stationary P1 z rho_mix xs

