% GenBif.m  Generates Bifurcation Diagram for Simplified slug model
%
% Written by Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
% Date 15.12.03
%
%
clear all
global data IX IU IFIXU uorg
%% Inputs for initialization
mG_in=0.362;
mL_in=8.64;
z=0.13;
P1_stasj=70.95;
h1_stasj=0.11;
P2_stasj=53.4;
P0=50; %Bar
%% Initialization
[x0,y_stasj,data]=initialize(z,mG_in,mL_in,P1_stasj,h1_stasj,P2_stasj,P0);

%% Defining simulation parameters and time vectors for open loop simulation
z_vec=[0.07;0.08;0.09;0.1;0.11;0.12;0.13;0.14;0.15;0.175;0.2;0.25;0.3;0.4;0.5;0.6;0.8;1];
Tvec=[12;12;12;24;48;70;70;15;10;5;3;3;3;3;3;3;3;3]*3600;
options=odeset('AbsTol',1e-6,'RelTol',1e-6,'MaxStep',60);

%% Simulation and Peak value identification
for i =1:length(z_vec)
    u=[z_vec(i),mG_in,mL_in,P0];simlength=Tvec(i);
    [t,x]=ode23t(@slug_3D_sim,[0;simlength],x0,options,u,'derivatives',data);
    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)];
    J=find(abs(t(I)-simlength*20/24)==min(abs((t(I)-simlength*20/24))));
    for j=J:length(I)
        y(j,:)=slug_3D_sim(t(I(j)),x(I(j),:),u,'measurements',data);
    end
    maxP1m(i)=max(y(J:length(t(I)),1));
    minP1m(i)=min(y(J:length(t(I)),1));
    maxP2m(i)=max(y(J:length(t(I)),2));
    minP2m(i)=min(y(J:length(t(I)),2));
    disp(sprintf('Simulation %i of %i completed',i,length(z_vec)))
end
z_vecm=z_vec;
save slugdata z_vecm maxP1m minP1m maxP2m minP2m

%% Setting Parameters for the tracking of the stationary solution
IX=1:3;
IU=4:7;
IFIXU=2:4;
%% Finding starting points for Broyden_secant
z=0.07;
u=[z;mG_in;mL_in;P0];
[t,x]=ode23t(@slug_3D_sim,[0;12*3600],x0,options,u,'derivatives',data);
z1=[x(end,:)';u];
z=0.075;
u=[z;mG_in;mL_in;P0];
[t,x]=ode23t(@slug_3D_sim,[0;12*3600],x0,options,u,'derivatives',data);
z2=[x(end,:)';u];uorg=u;
%% Tracking stationary solution
[out,uu,step1,tt]=secant_cont_circle('slug_wrapper',z2,z1,0.05);
%% Computing measurements for stationary solution
disp('Computing measurements')
xs=out(:,IX)';us=out(:,IU)';
for i=1:length(out)
    ys(i,:)=slug_3D_sim(0,xs(:,i),us(:,i),'measurements',data);
end
P1=ys(:,1);
P2=ys(:,2);
z=us(1,:);
save Stationary P1 P2 z xs

%% Plotting Bifurcation chart together with OLGA data
if exist([pwd '\OLGAslugdata.mat'])==2
    load OLGAslugdata
    figure(1)
    plot(us(1,:),ys(:,1),'r--',z_vecm,maxP1m,'r-',z_vecm,minP1m,'r-',z_vec,maxP1/1e5,'b-',z_vec,minP1/1e5,'b-',z_stasj,P1stasj,'bx')

    figure(2)
    plot(us(1,:),ys(:,2),'r--',z_vecm,maxP2m,'r-',z_vecm,minP2m,'r-',z_vec,maxP2/1e5-50,'b-',z_vec,minP2/1e5-50,'b-',z_stasj,P2stasj-50,'bx')

else %plotting only Matlab data
    figure(1)
    plot(us(1,:),ys(:,1),'r--',z_vecm,maxP1m,'r-',z_vecm,minP1m)

    figure(2)
    plot(us(1,:),ys(:,2),'r--',z_vecm,maxP2m,'r-',z_vecm,minP2m)
end
    
    