clear all
mG_in=6.1717e-5*1.5;
P_feed=0.31e5+1.013e5;
z=0.2;uLs=0.1528;
P1_stasj=0.127e5+1.013e5;
h1_stasj=0.019/2*0.6;%H1=0.2128
P2_stasj=0.037e5+1.013e5;
P0=1.013e5;
gen_data
%Initialization
[x0,y_stasj,data]=initialize(z,mG_in,P_feed,P1_stasj,h1_stasj,P2_stasj,P0,uLs);


z_vec=[0.05;0.08;0.12;0.14;0.16;0.18;0.19;0.20;0.21;0.22;0.24;0.26;0.3;0.35;0.4;0.5;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,P_feed,P0];simlength=600;
    [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 slugdata2 z_vecm maxP1m minP1m maxP2m minP2m
