clear all
mG_in=0.362;
mL_in=8.64;
z=0.13;
P1_stasj=71.7e5;
h1_stasj=0.1;
P2_stasj=53.5e5;
%Initialization
[x0,y_stasj,data]=initialize(z,mG_in,mL_in,P1_stasj,h1_stasj,P2_stasj);
z_vec=[0.05;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];
z_vec=1;
M=eye(4);M(4,4)=0;
options=odeset('Mass',M,'BDF','on','AbsTol',1e-8,'RelTol',1e-8,'MaxStep',60);
x0=[x0;data.rho_mix_i];
simlength=10*3600;

for i =1:length(z_vec)
    u=[z_vec(i),mG_in,mL_in];
    [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*2/3)==min(abs((t(I)-simlength*2/3))))
    maxP1(i)=max(y(J:length(t(I)),1))
    minP1(i)=min(y(J:length(t(I)),1))
    maxP2(i)=max(y(J:length(t(I)),2))
    minP2(i)=min(y(J:length(t(I)),2))
    x0=x(end,:);
    
    i
end
save slugdata z_vec maxP1 minP1 maxP2 minP2