% Simulation of the new 4d model
clc
clear all
% **** System parameters ****
par = v4_new_4d_parameters();
% **** System Inputs ****
z0 = 0.5;
%Qg_in = 0.0023569; %[m3/sec]
wG_in0 = 2e-3;%0.0023569;%Qg_in*par.Rog_n; %[kg/sec]
QL_in = 0.4;%0.39202239; %[L/sec]
wL_in0 = QL_in*par.rho_w/1000; %[kg/sec]
global u0;
u0 = [z0;wG_in0;wL_in0];
% ******** Finding initial states *************
disp('Initializing...')
tic
[x0,y0,par] =  v4_new_4d_initialize(u0,par);
toc
disp('Initialization completed.');
disp('#########################');

%dt = 0.005;    % Simulation step time (sec)
t_f = 2;        % Final simulation time (hour)

%% **********************
disp('Simulation of nonlinear model');
dx = [-0.01;0;0;0];
u=[z0;wG_in0;wL_in0];
M=eye(4);
options=odeset('Mass',M,'BDF','on','AbsTol',1e-8,'RelTol',1e-8,'MaxStep',1);
[t,x]=ode15s(@v4_new_4d_model,[0;3600*t_f],x0 + dx,options,u0,'derivatives',par);
% *** Calculating measurements ****
n = length(t);
y = zeros(3,n);
for k = 1:n
    yt = v4_new_4d_model(t(k),x(k,:),u0,'measurements',par);
    y(:,k) = yt([1:2 4]);
end

disp('Simulation completed.');
disp('#########################');

figure(2)
clf
set(gcf,'Color',[1,1,1])
set(gca,'FontSize',8)

subplot(3,1,1)
set(gca,'FontSize',8,'OuterPosition',[0  2/3  1  1/3])
plot(t,y(1,:),'k');
title('a) Pressure at inlet of pipeline (\itP\rm\bf_1)','FontSize',10,'FontWeight','bold','FontName','Times');
ylabel('\itP\rm_1 [bar]','FontSize',10,'FontName','Times')
xlabel('time [sec]','FontSize',10,'FontName','Times')
xlim([0 t_f*3600])
%ylim([68 80])
min(y(1,:))
max(y(1,:))

subplot(3,1,2)
set(gca,'FontSize',8,'OuterPosition',[0  1/3  1  1/3])
plot(t,y(2,:),'k');
title('b) Pressure at top of riser (\itP\rm\bf_2)','FontSize',10,'FontWeight','bold','FontName','Times');
ylabel('\itP\rm_2 [bar]','FontSize',10,'FontName','Times')
xlabel('time [sec]','FontSize',10,'FontName','Times')
xlim([0 t_f*3600])
%ylim([50 62])


subplot(3,1,3)
set(gca,'FontSize',8,'OuterPosition',[0  0  1  1/3])
plot(t,y(3,:),'k');
title('c) Mass flow rate form choke valve (\itW_{out}\rm\bf)','FontSize',10,'FontWeight','bold','FontName','Times');
ylabel('\itW_{out}\rm [kg/s]','FontSize',10,'FontName','Times')
xlabel('time [sec]','FontSize',10,'FontName','Times')
xlim([0 t_f*3600])
%ylim([5 35])
