% GenBif.m  Generates Bifurcation Diagram for Simplified slug model
%
% Written by Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
% Date 15.12.03
%
%
%clear all
clear
global data IX IU IFIXU uorg
%% Inputs for initialization
Exp_data;


%% Initialization
[x0,y_stasj,data]=initialize(z,mG_in,mL_in,P1_stasj,h1_stasj,P2_stasj);

%% Defining simulation parameters and time vectors for open loop simulation
z_vec=[0.07;0.1;0.13;0.14;0.15;0.16;0.17;0.18;0.20;0.22;0.25;0.28;0.3;0.35;0.5;0.6;0.7;1];
Tvec=[12;12;12;24;48;48;48;48;48;48;48;48;48;24;24;24;24;24]*10;
options=odeset('AbsTol',1e-8,'RelTol',1e-8,'MaxStep',60);

%% Simulation and Peak value identification
for i =1:length(z_vec)
    u=[z_vec(i);mG_in;mL_in]; 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))));
    %clear y %purpose of this ????
    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));
    %maxP1m(i)=max(y(J:end,1));
    %minP1m(i)=min(y(J:end,1));
    %maxP2m(i)=max(y(J:end,2));
    %minP2m(i)=min(y(J:end,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:6;
IFIXU=2:3;
%% Finding starting points for Broyden_secant
z=0.1;
u=[z;mG_in;mL_in];
[t,x]=ode23t(@slug_3D_sim,[0;12*600],x0,options,u,'derivatives',data);
z1=[x(end,:)';u];
z=0.13;
u=[z;mG_in;mL_in];
[t,x]=ode23t(@slug_3D_sim,[0;12*600],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.001);
%% 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)-1.013;
P2=ys(:,2)-1.013;
z=us(1,:);
save Stationary P1 P2 z xs
save bifurksim z_vecm maxP1m minP1m maxP2m minP2m

%% Plotting Bifurcation chart together with OLGA data
if exist([pwd '\eksperimentelle_data.mat'])==2
    load Stationary
    load eksperimentelle_data
    load bifurksim
    subplot(2,1,1)
    %plot(z,P1,'r--',z_vecm,maxP1m-1.013,'r-',z_vecm,minP1m-1.013,'r-',Z/100,P1_max,'b:',Z/100,P1_min,'b:')
    %legend('Stationary solution','Simulated values','','Experimental values','')
    plot(z_vecm,maxP1m-1.013,'r-',Z/100,P1_min,'b*',z_vecm,minP1m-1.013,'r-',Z/100,P1_max,'b*',z,P1,'r--')
    axis([0 1 0 0.6])
    ylabel('P_1 [bara]')
    legend('Simulated values', 'Experimental values')
    subplot(2,1,2)
    plot(z,P2,'r--',z_vecm,maxP2m-1.013,'r-',z_vecm,minP2m-1.013,'r-',Z/100,P2_max,'b*',Z/100,P2_min,'b*')
    %legend('Stationary solution','Simulated values','','Experimental values','')
    axis([0 1 0 0.4])
    ylabel('P_2 [bara]')
else %plotting only Matlab data
    figure(1)
    grid on
    plot(z(1,:),P1,'b--',z_vecm,maxP1m,'r-',z_vecm,minP1m,'r-')

    figure(2)
    plot(z(1,:),P2,'b--',z_vecm,maxP2m,'r-',z_vecm,minP2m,'r-')
end
    
    