clear all
mG_in=0.362;
mL_in=8.64;
z=0.13;
P1_stasj=70.13e5;
h1_stasj=0.1;
P2_stasj=53.4e5;
P0=50e5;
[x0,y_stasj,data]=initialize(z,mG_in,mL_in,P1_stasj,h1_stasj,P2_stasj,P0);
z_vec=[0.07;0.08;0.09;0.1;0.11;0.12;0.13;0.14;0.15;0.175;0.2;0.25;0.3;0.4;0.5;0.6;0.8;1];
Tvec=[12*3600;12*3600;12*3600;24*3600;48*3600;70*3600;70*3600;15*3600;10*3600;5*3600;3*3600;3*3600;3*3600;3*3600;3*3600;3*3600;3*3600;3*3600];
M=eye(4);M(4,4)=0;
options=odeset('Mass',M,'BDF','on','AbsTol',1e-8,'RelTol',1e-8,'MaxStep',100);

x0=[x0;data.rho_mix_i];
for i =1:length(z_vec)
    u=[z_vec(i),mG_in,mL_in,P0];simlength=Tvec(i);
    [t,x]=ode23t(@slug_3DAE,[0;simlength],x0,options,u,'derivatives',data);
    I=[1:10:10*((rem(length(t),10)~=0)*(floor(length(t)/10)+1)+(rem(length(t),10)==0)*(floor(length(t)/10))),length(t)];
    for j=1:length(I)
        y(j,:)=slug_3DAE(t(I(j)),x(I(j),:),u,'measurements',data);
    end
    J=find(abs(t(I)-simlength*20/24)==min(abs((t(I)-simlength*20/24))))
    maxP1m(i)=max(y(J:length(t(I)),1))
    minP1m(i)=min(y(J:length(t(I)),1))
    maxP2m(i)=max(y(J:length(t(I)),2))
    minP2m(i)=min(y(J:length(t(I)),2))
    x0=x(end,:);
    i
end
z_vecm=z_vec;
save slugdata z_vecm maxP1m minP1m maxP2m minP2m
Track_Stationary