%% Riser-Annulus Model Simulation
clear all;
clc;
global u0;
%% Inputs to model 
% Manipulated variables
z1=0.6;    % production choke opening (u1) 
z2=0.4;    % injection valve opening(u2) 
% Disturbances
GOR = 0;                     % Mass GOR from reservoir, NB: MASS, increasing gives: lower BH, stable
P_r = 160;                      % Rservoir pressure [bar]
P_gs = 140;                     % Gass source pressure [bar]
P0 = 15;                        % Separator Pressure [bar]
d = [GOR;P_r;P_gs;P0];         % Disturbance Vector
u0=[z1;z2;d]; %Only u(1) and u(2) are MVs

%% ***** Parameter Constants*******
par = v1_gaslift_parameters();

%% Initial States 
disp('************ Start search for equilibrium point *****************');
tic
[x0,flag] = v1_gaslift_initialize(u0,par);
toc
disp('************ End search for equilibrium point *****************');
%% Linearazation linmod, Matlab linerazation on simulink model
disp('************ Linearization *****************');
    %[A B C D]=central_difference(x0,u0,par); % State Space Model
    [A,B,C,D]=linmod('linerazation_simulink',x0,u0);
    [eigval]=real(eig(A));
    unstable=0;
    for i=1:length(eigval)
        if eigval(i) >=0
            unstable=1;
        end
    end
    if(unstable)
        disp('UNSTABLE')
    else
        disp('STABLE')
    end
dx=[-1;0;0;];

%% Simulation, solve diff.equation 
disp('********** Start Simulation of nonlinear model *******************');
M=eye(3);
T_span = 1:10*3600;       
%options=odeset('Mass',M,'AbsTol',1e-12,'RelTol',1e-12,'MaxStep',1,'BDF','on','Maxorder',2);
options=odeset('AbsTol',1e-10,'RelTol',1e-10,'MaxStep',100,'Mass',M);
[t,x]=ode15s(@v1_gaslift_model,T_span,x0+dx,options,u0,'derivatives',par);
disp('********** End Simulation of nonlinear model *******************');

%% Measurements 
n = length(t);
for k = 1:n
    y_k = v1_gaslift_model(t(k),x(k,:),u0,'measurements',par);
    y_k=y_k';
    y(k,:) = y_k;
end

%% Plotting
plotting(x,y,t);