%script for aa beregne xinit(=x0) for ulike ventilaapninger for bruk i track stationary.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc
% Inputs for initialization
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);

%Inputs for simulation
z=0.16;
u=[z;mG_in;P_feed;P0]
x0=[x0;data.rho_mix_i];
%Setting integrator options
M=eye(4);M(4,4)=0;
options=odeset('Mass',M,'BDF','on','AbsTol',1e-8,'RelTol',1e-8,'MaxStep',60);
%Simulating
disp('Simulating')
[t,x]=ode23t(@slug_3DAE,[0;1000],x0,options,u,'derivatives',data);
disp('Computing measurements')
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 i=1:length(I)
    y(i,:)=slug_3DAE(t(I(i)),x(I(i),:),u,'measurements',data);
end
t=t/60;
figure(1)
title('Simplified slug model')
plot(t(I),y(:,1))
xlabel('time [min]')
ylabel('Upstream pressure P_1 [Bar]')
x0=x(end,:);
save xinit016 x0
save data data