% Simplified model of twophse floe in pipeline with slug behaviour
% 16.05.04
% Gisle Otto Eikrem
% Inkludert forskjellige inflow betingelser av Henrik Tuvnes 19.12.08

close all;
clear all;

% Flow from reservoir options
%flow_type =  'constant ';         % Constant inflow
flow_type = 'linear   ';           % Linear inflow, w=PI*(Pres - Pbh)%
%flow_type = 'reservoir';        % Dynamic reservoir inflow

GOR_type='constant';
%GOR_type='linear';

controller_type = 'openloop'; % Open loop simulation
%controller_type = 'dhpcontr'; % Closed loop simulation

%format short e;

Ts_controller = 1; % Time step for execution of controller
min = 1000;
Ts = 1;
Num = Ts_controller/Ts; %Number of iterations between each execution of the controller
L = min*60/Ts;

u_PC_simple(1:Num-1) = 0.4; % The initial process input sequence
p_bottom_riser_bara(1) = 4;

C0 = 400*60/Ts; %Choking back at this timestep
C1 = 800*60/Ts; % Closing the loop at this timestep
C2 = 1500*60/Ts; %Opening the loop at this timestep


matrix_init; % Initialize the sizes of the matrices

%load pipeline_model.mat;

eikrem_pipeline_geometry; % load the geometry of the selected pipeline

eikrem_init; % Initial values of the simplified pipeline model

if flow_type == 'reservoir',
reservoir_init;   % initialize reservoir parameters and conditions
else end

for i=1:L+Num-1, % Iteration loop
 
 time_log(i) = i*Ts; % Log of the time vector
 
    if mod(i,Num)==0, % The controller only executed at the given timestep
       controller_pipeline;
    end     % Execution of controller
 
% Calculate the change of mass for the next time step, i+1
  eikrem_pipeline_model; 
  m(:,i+1) = m(:,i) + Ts*dm; % calculate the a priori estimate of the state vector 
  
  if m(1,i+1) < 0.001, % This inequality prevents the estimated mass to drop below the 
      m(1,i+1) = 0.001;
      error_1 = 1
  end
 
  if m(2,i+1) < 0.001, % This inequality prevents the estimated mass to drop below the 
      m(2,i+1) = 0.001;
      error_2 = 1      
  end
  
  if m(3,i+1) < 0.001, % This inequality prevents the estimated mass to drop below the 
      m(3,i+1) = 0.001;
      error_3 = 1      
  end
  
  if m(4,i+1) < 0.001, % This inequality prevents the estimated mass to drop below the 
      m(4,i+1) = 0.001;
      error_4 = 1
  end

  if mod(i,100)==0, disp([num2str(floor(i/L*100)) '/100%']), end
  
end % Iteration loop

m=m(1:end,1:end-1); % remove last calc

hours=time_log./3600;  % simulated time in hours

save('pipeline_model.mat');

plot_pipeline;
%plot_report;