%% prod_vs_injectionrate_direct_import
clear all; clc;
dW=28; % Number of different gas injections rates pr valve openings 
dS=10; % Number of different gas injections rates pr valve openings in stability Map
par=v1_gaslift_parameters();

% Inputs to model 
% Manipulated variables
z1=0.75;    % production choke opening (u1) 
z2=1.25;    % 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 = 150;                     % 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

%% Generate Data Vector
simData=[];
stabData=[];
for z1=0.1:0.1:1 % u(1)
    for w_g=0.3:0.1:3
        disp([' z1: ' , num2str(z1) , ' w_G: ', num2str(w_g) ])        
        u0(1)=z1;
        u0(2)=w_g;
        [x_eq,flag,y] = v1_gaslift_initialize(u0,par);
        prod_valve_opening=z1;
        gas_injection_rate=w_g;
        production_rate=y(6);    
        simData=[simData; prod_valve_opening, gas_injection_rate, production_rate];
        
        stable=1;
        [A B C D]=central_difference(x_eq,u0,par); % State Space Model
        [eigval]=real(eig(A));
        for i=1:length(eigval)
            if eigval(i) >=0
                stable=0;
            end
        end            
        stabData=[stabData; u0(1) u0(2) stable] ;
    end
end
% Save Generated Data
save('matlab_prod_vs_gasinjecton.mat','simData')
save('matlab_stab_data.mat','stabData');

%% Load Old Data Instead of Generating New
var=load('matlab_prod_vs_gasinjecton.mat');
simData=var.simData;
var2=load('matlab_stab_data.mat');
stabData=var2.stabData;

%% Import StabilityMap made by different file
stabTrans=[];%find the stability transition points
for i=1:dW:size(stabData)
    
    vec=stabData(i:i+dW-1,:);
    
    for j=1:size(vec,1)-1
        if vec(j,3)==0 && vec(j+1,3)==1
            stabTrans=[stabTrans; vec(j+1,1:2)];
            break;
        end
    end
end
stabTrans=stabTrans(:,2);


%% Plotting Production vs gas injection
a=size(simData);
str0 =['- -'];
str2= ['\leftarrow'];
teller=1;

for i=1:dW:size(simData)    
    fig=figure(1);
    x=simData(i:i+dW-1,2);
    y=simData(i:i+dW-1,3);
    
    plot(x,y,'k')
    hold on;
    plot(x,y,'+k')
    hold on;
    
    for k=1:1:length(x)
        if isequal(num2str(x(k)),num2str(stabTrans(teller)))
        %if x(k) == stabTrans(teller)            
            trans=plot(stabTrans(teller),y(k),'rx','MarkerSize',15,'LineWidth',2);
            hold on;
            break;
        end
    end
    teller=teller+1;
 
        if i>=dW*7 && i<=dW*9
            for h=1:1:power(2,(a-i)/dW)                   
                str2= [str2, str0];                
            end
        end
        
        str1 = [str2 ' u = ' num2str(100*simData(i,1)) '%'];
        txt=text(simData(i+dW-1,2),simData(i+dW-1,3),str1,'HorizontalAlignment','left');
        set(txt,'Color','red','FontSize',15);
        title('Low-Order Model Steady-State Production Flow Rate  vs Gas Injection Rate','FontSize',20,'FontWeight','bold')
        ylabel('Production Flow [kg/s]','FontSize',20,'FontWeight','bold')
        xlabel('Injected Gas to Annulus [kg/s]','FontSize',20,'FontWeight','bold')
        set(gca,'FontSize',16);
        %ylim([11 18])
        xlim([simData(1,2) 1.25*simData(dW,2)])   
        
    str2= ['\leftarrow'];
end

set(gca,'xTick',0:0.2:4)
set(gca,'yTick',4:1:28)
legend(trans,'Transition To Stability');

