%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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
Exp_data;
%Initialization
[x0,y_stasj,data]=initialize2(z,mG_in,mL_in,P1_stasj,h1_stasj,P2_stasj);


%Inputs for simulation
z=0.30;
u=[z;mG_in;mL_in];
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;3600*6],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_3DAE(t(I(i)),x(I(i),:),u,'measurements',data);
end

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