%Track_stationary.m
clear all
global data IX IU IFIXU uorg
load data
IX=1:4;
IU=5:8;
IFIXU=2:4;
mG_in=0.362;
mL_in=8.64;
P0=50e5;
z=0.07;
load xinit007
u=[z;mG_in;mL_in;P0];
z1=[x00;u];
load xinit0075
z=0.075
u=[z;mG_in;mL_in;P0];
uorg=u;
z2=[x0;u];
[out,uu,step1,tt]=secant_cont_circle('slug_wrapper',z2,z1,0.05);
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 OLGAslugdata

figure(1)
plot(us(1,:),y(:,1),'r--',z_vecm,maxP1m,'r-',z_vecm,minP1m,'r-',z_vec,maxPrb/1e5,'b-',z_vec,minPrb/1e5,'b-',z_stasj,Prbstasj,'bx')
axis([0 1 60 80])
figure(2)
plot(us(1,:),y(:,2),'r--',z_vecm,maxP2m,'r-',z_vecm,minP2m,'r-',z_vec,maxP2/1e5,'b-',z_vec,minP2/1e5,'b-',z_stasj,P2stasj,'bx')
axis([0 1 50 60])
P1=y(:,1);
P2=y(:,2);
z=us(1,:);
h1=xs(4,:);
save Stationary P1 P2 z h1 xs
% plot(us(1,:),y(:,3))
% P1=y(:,3);
% z=us(1,:);
% rho_mix=xs(4,:);rho_T=xs(5,:);
% save Stationary4 P1 z rho_mix rho_T xs
