% 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
uL=2.39; %vannmengde inn l/min
uG=1.58; %luftmengde inn l/min
P_feed=1.325;
z=0.19;
uLs=uL/1000/60/pi/(0.0195/2)^2;
P1_stasj=1.133;
h1_stasj=0.019/2*0.6;
P2_stasj=1.036;
P0=1.013; %Bar
temp=293;


Ml=28.84/1000; %Molarmasse luft [kg/mol]
uLs=uL/1000/60/pi/(0.0195/2)^2; %Water rate (m/S)  %må endres
R=0.08205783*1.013; %l bar K^-1 mol^-1
rho_in= P1_stasj*Ml/R/temp;   %[kg/l]
mG_in= rho_in*uG/60; %(6.1717e-5)*1.5; % [kg/s] Endres 


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

%% Defining simulation parameters and time vectors for open loop simulation
z_vec=[0.07;0.1;0.13;0.16;0.17;0.18;0.21;0.22;0.23;0.24;0.25;0.28;0.3;0.4;0.5;0.6;0.8;1];
Tvec=[12;12;12;24;48;48;48;48;48;20;20;20;20;12;3;3;3;3]*100;
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,P_feed,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))));
    %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:7;
IFIXU=2:4;
%% Finding starting points for Broyden_secant
z=0.1;
u=[z;mG_in;P_feed;P0];
[t,x]=ode23t(@slug_3D_sim,[0;12*600],x0,options,u,'derivatives',data);
z1=[x(end,:)';u];
z=0.11;
u=[z;mG_in;P_feed;P0];
[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);
P2=ys(:,2);
z=us(1,:);
save Stationary P1 P2 z xs

%% Plotting Bifurcation chart together with OLGA data
if exist([pwd '\eksperimentelle_data.mat'])==2
    load eksperimentelle_data
    figure(1)
    plot(us(1,:),ys(:,1),'r--',z_vecm,maxP1m,'r-',z_vecm,minP1m,'r-',Z/100,P1_max+1.013,'b.',Z/100,P1_min+1.013,'b.')

    figure(2)
    plot(us(1,:),ys(:,2),'r--',z_vecm,maxP2m,'r-',z_vecm,minP2m,'r-',Z/100,P2_max,'b-',Z/100,P2_min,'b-')

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

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