%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   slug_3D_script.m
%
%   Script for simulating the simplified slug model in open loop
%
%   Written by :    Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
%   Date:           08.10.02
%   
%   For questions, comments or documentation, contact author
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
%% Inputs for initialization
mG_in=6.17e-5*1.5;
P_feed=1.325;
z=0.3;
uLs=2.37/1000/60/pi/(0.0195/2)^2;
P1_stasj=1.145;
h1_stasj=0.019/2*0.6;
P2_stasj=1.038;
P0=1.013; %Bar
%% Initialization
[x0,y_stasj,data]=initialize(z,mG_in,P1_stasj,h1_stasj,P2_stasj,P0,P_feed,uLs);


%% Simulation Boundary data
z=0.29;
u=[z;mG_in;P_feed;P0];
x0=[x0];

options=odeset('AbsTol',1e-8,'RelTol',1e-8,'MaxStep',60);
%% Simulating
disp('Simulating')
[t,x]=ode23t(@slug_3D_sim,[0;600],x0,options,u,'derivatives',data);
%% Computing Measurements
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_3D_sim(t(I(i)),x(I(i),:),u,'measurements',data);
end

%% Displaying results
t=t/60;
figure(1)
title('Simplified slug model')
plot(t(I),y(:,1))
xlabel('time [min]')
ylabel('Upstream pressure P_1 [Bar]')
figure(2)
plot(t(I),y(:,2));
title('Simplified slug model')
xlabel('time [min]')
ylabel('Downstream pressure P_2 [Bar]')