function[initComb,     ConAnComb,...
         initCombNN,   ConAnCombNN]...
          = InitOpt(par,scl,con,ther,def,itg,dirmain)

% This function returns initial values and annonymous functions for the 
% submodels.

% Written by Julian Straus, Trondheim 09.12.2018

% The input structures are defined as
%       par:    Parameters used in the calculation, output of the system
%       scl:    Scaling factors of the system, must include:
%               .Q      Scaling factor for the flow
%               .p      Scaling factor for the pressure
%               .T      Scaling facto for the temperature
%               .W      Scaling factor for the pwer
%       con:    Physical constants, must contain
%               .R      Ideal Gas constant
%       ther:   Thermodynamical values to calculate henry coefficients of
%               the supercritical components, must include:
%               .Hi    	Henry coefficient, i = 1-3
%       def:    Defined stream values, e.g pressure and composition for all
%               inlets and the pressure of the outlets or the flow of the
%               inlets
%       itg:    Intiger values, which should include
%               .k    	Total number of variables per stream
%               .j     	Total number of components in each stream


% Load of the initial values from the flowsheet
load([dirmain,'\InitTot_FO_MakeUpSec.mat'])         % Solution for the MakeUp Section
load([dirmain,'\InitTot_FO_ReacSec_n_12=10.mat'])   % Solution for the Reaction Section
load([dirmain,'\InitTot_FO_SepSec.mat'])            % Solution for the Separation Section

% Addition of the path of the reduced order unit operations
UO = [dirmain,'\Util4. Unit Operations\Simulation'];
addpath(UO)

%% Reassignment of the variables

% Calculation of the values for the MakeUp section
x.Comp1 = [initMakeUpSaved(1:8); par.M{1}.S; 6.75e6; initMakeUpSaved(9:16)];
x.HEx1  = initMakeUpSaved(17:40);
x.Comp2 = [5.2307e6; initMakeUpSaved(41:48)];
x.Mix1  = initMakeUpSaved(49:64);
x.Comp3 = [3.0463e6; initMakeUpSaved(65:72)];

% Calculation of the values for the Reaction section
x.HEx2   = initReactSaved(1:32);
x.Split1 = [initReactSaved(33:56); def.R.Sratio{1}(1:2)];
x.Val1   = initReactSaved(57:65);
x.HEx3   = initReactSaved(66:90);
x.HEx4   = initReactSaved(91:114);
x.Mix2   = initReactSaved(115:122);
x.Reac1  = initReactSaved(123:194);
x.Reac2  = initReactSaved(195:266);
x.HEx5   = [def.R.Tref; initReactSaved(267:282); 3.7315e+07];

% Calculation of the values for the Separation section
x.HEx6   = initSepSaved(1:32);
x.Split2 = [initSepSaved(33:48); def.S.Sratio{1}(2)];
x.HEx7   = initSepSaved(49:73);
x.Split3 = [initSepSaved(74:89); def.S.Sratio{2}(2)];
x.HEx8   = [def.S.Tref(1); initSepSaved(90:113); 5.1456e6];
x.HEx9   = [def.S.Tref(2); initSepSaved(114:137); 4.4178e6];
x.Mix3   = initSepSaved(138:145);
x.HEx10  = [def.S.Tref(3); initSepSaved(146:169); 8.1642e6];
x.Sep    = initSepSaved(170:177);

%% Definition of the initial values for the different sections using the
% full order models

initMakeUp_In  = [x.Comp1; x.HEx1; x.Comp2; x.Mix1; x.Comp3];
initReact_In   = [x.HEx2; x.Split1; x.Val1; x.HEx3; x.HEx4; x.Mix2; x.Reac1; x.Reac2; x.HEx5];
initReactNN_In = [x.HEx2(9:24); x.Split1(end-1:end); x.HEx5([1 end])];
initSep_In     = [x.HEx6; x.Split2; x.HEx7; x.Split3; x.HEx8; x.HEx9; x.Mix3; x.HEx10; x.Sep];
initSepNN_In   = [x.HEx6([1:8 17:24]); x.Sep; x.HEx7(17:24); x.HEx6(25:32);...
                  def.S.Sratio{1}(2); def.S.Tref';... Parameters
                  x.HEx8(end); x.HEx9(end); x.HEx10(end);  x.Split3];
initComb_In    = [initMakeUp_In; initReact_In; initSep_In];

% Rescaling of the initial values
initMakeUp  = initMakeUp_In./scl.M{11}';
initReactNN = initReactNN_In./scl.R{23}';
initSepNN   = initSepNN_In./scl.S{23}';
initComb    = initComb_In./[scl.M{11} scl.R{19} scl.S{19}]';
initCombNN  = [initMakeUp; initReactNN; initSepNN];
        
        
%% Definition of the anonymous functions for the combined flowsheet

% Original model
itg.nvar    = numel(initComb);
def.M.fun      = @(nvar,into,colind) DefMakeUp_Comb(nvar,into,colind,scl.M,def.M);
def.R.fun      = @(nvar,into,colind) DefReact_Comb(nvar,into,colind,scl.R,def.R);
def.S.fun      = @(nvar,into,colind) DefSep_Comb(nvar,into,colind,scl.S,def.S);
def.M.ineqfun  = @(nvar,into,colind) InEqMakeUp(nvar,into,colind,scl.M,def.M);
% def.R.ineqfun  = @(nvar,into,colind) InEqReact(nvar,into,colind,scl.R,def.R);
% def.S.ineqfun  = @(nvar,into,colind) InEqSep(nvar,into,colind,scl.S,def.S);
def.R.ineqfun  = @(nvar,into,colind) InEqReact_Un(nvar,into,colind,scl.R,def.R);
def.S.ineqfun  = @(nvar,into,colind) InEqSep_Un(nvar,into,colind,scl.S,def.S);
ConAnComb      = @(x) ConComb_min(x,par,scl,con,ther,def,itg);

% Neural network
itg.nvar    = numel(initCombNN);
def.R.fun       = @(nvar,into,colind) DefReact_NN_Comb(nvar,into,colind,scl.R,def.R);
def.S.fun       = @(nvar,into,colind) DefSep_NN_Comb(nvar,into,colind,scl.S,def.S);
% def.R.ineqfun   = @(nvar,into,colind) InEqReact_NN(nvar,into,colind,scl.R,def.R);
% def.S.ineqfun   = @(nvar,into,colind) InEqSep_NN(nvar,into,colind,scl.S,def.S);
def.S.ineqfun   = @(nvar,into,colind) InEqSep_NN_Un(nvar,into,colind,scl.S,def.S);
def.R.ineqfun   = @(nvar,into,colind) InEqReact_NN_Un(nvar,into,colind,scl.R,def.R);
ConAnCombNN     = @(x) ConComb_NN_min(x,par,scl,con,ther,def,itg);
