clc
clear all
% **** System Inputs ****
par = v4_new_4d_parameters();

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]

%QG_in = 4.5; %litr/min
%QL_lin = 4; %litr/min
%wG_in0 = par.Rog_n*QG_in/(60*1000); %kg/s
%wL_in0 = par.rho_L*QL_lin/(60*1000); %kg/s
% **** Automatic ********
 zs = [0.2:0.02:0.24 0.26:0.01:0.30 0.32 0.35 0.40:0.1:1];
ns0 = size(zs);
ns = ns0(2);
y1min = zeros(ns,1);
y1max = zeros(ns,1);
y1init = zeros(ns,1);
y1end = zeros(ns,1);

y2min = zeros(ns,1);
y2max = zeros(ns,1);
y2init = zeros(ns,1);
y2end = zeros(ns,1);

y3min = zeros(ns,1);
y3max = zeros(ns,1);
y3init = zeros(ns,1);
y3end = zeros(ns,1);

dx = [-0.0001;0;0;0];
M=eye(4);
options=odeset('Mass',M,'BDF','on','AbsTol',1e-12,'RelTol',1e-12,'MaxStep',1);

for s = 1:ns
% **** System Inputs ****
z0 = zs(s);
disp(['Z=' num2str(z0)])
u0 = [z0;wG_in0;wL_in0];
% ******** Initializing *************
[x0,y0,par] = v4_new_4d_initialize(u0,par);
% **********************
if (z0<0.12)
    tf = 4;      
else
    tf = 4;  
end

u=[zs(s);wG_in0;wL_in0];

%** Simulation using Simple method *****
[t,x]=ode15s(@v4_new_4d_model,[0;3600*tf],x0 + dx,options,u,'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:3);
    end
%n1 = floor(n/2);
y1min(s) = min(y(1,1:n));%n1:n));
y1max(s) = max(y(1,1:n));%n1:n));
y1init(s) = y(1,1);
y1end(s) = y(1,end);

y2min(s) = min(y(2,1:n));%n1:n));
y2max(s) = max(y(2,1:n));%n1:n));
y2init(s) = y(2,1);
y2end(s) = y(2,end);

y3min(s) = min(y(3,1:n));%n1:n));
y3max(s) = max(y(3,1:n));%n1:n));
y3init(s) = y(3,1);
y3end(s) = y(3,end);
clear y x t;
end

save('matlab_bif_data.mat','y1min','y1max','y1init','y1end','y2min','y2max','y2init','y2end','y3min','y3max','y3init','y3end');


Bif_Plot2;
%Steady_State_auto;