%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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.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.30;
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;200],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


subplot(2,2,1)
plot(t(I),y(:,1)-1.013)
ylabel('P_1 [barg]')

subplot(2,2,2)
plot(t(I),y(:,2)-1.013)
ylabel('P_2 [barg]')

subplot(2,2,3)
plot(t(I),y(:,4))
ylabel('W [kg/s]')
xlabel('tid [s]')

subplot(2,2,4)
plot(t(I),y(:,5))
ylabel('Q [10^{-4} m^3/s]')
xlabel('tid [s]')
% g=gcf;
% set(g,'PaperUnits','centimeters','PaperType','A4')
% opts = struct('FontMode','fixed','FontSize',9,'Width',15,'Color','rgb');
% exportfig(g,'simulering_z_40.eps',opts)




%  sys=[P1/1e5,P2/1e5,rho_T,mdot_mix,mdot_mix/rho_T];