% This is a script file to calculate the second derivative of cost with
% respect to the unconstrained degrees of freedom Reflux , Side stream and 
% the liquid split ratio(L, S and Rl)
%%Three Product Petlyuk Column
%Constant Relative Voltality
% Initializes the model and runs
%%
clear all; close all; clc;
global coldata dyndata hydraulic termodata UU alfa tend NI i
%% COLUMN DATA %%
% Number of stages in 6 sections of Petlyuk 
coldata.Nt = [12 12 12 12 12 12];
% Initial holdup stages
coldata.M0 = [1 1 1 1 1 1 1 1 1 1 1 1];
%Boiling points of ethanol, propanol and n-butanol
coldata.Tboil = [78.37 97 117.4];
termodata.P = 1.0;
hydraulic.k1 = 1/0.063;
%% THERMODYNAMIC DATA %%
alfa=[4.215374389	2.102194688	1] %ethanol propanol butanol
%delta=1e-5;

%% load nominal x u
load Nominal_final.mat        % The nominal value obtained from optimization for ...
%minimizing the sum of impurities in the product streams (D, S and R)
UU=[];                        %
i=0;%to initiate Petlyuk_openloop make change in CS by i+1
indmat_Yes_S1
% make nominal x and u global
dyndata.x_nom = x_nom; 
dyndata.u_nom=u_nom;%           we have made initial input global
dyndata.temp_nom=temp_nom;
%% Run this file to get initial values with LV Configuartion, uncomment
% line number 
NI=5; % Number of integerators
options = odeset('MaxStep',10);
tend=20000;
[t,x]=ode15s(@Simulation_Petlyuk,[0 tend],[dyndata.x_nom;zeros(NI,1)],options);
%% xf states to plot profiles.
xf=x(end,1:225)';
u_nom=UU(end,2:12)';
%dxdt = norm(Petlyuk_dxdt(0, xf, u_nom));
% Print the values of compositions in D (distillate), S (side stream), R (residue)...
%and D1 and R1, (top and bottom product from the prefractionator both  end)
%dis_composition
%result_simulation = sprintf (strcat(' min sum of impurities\n',...
  %  'distillate flow rate,         D[kmol/min]      = %1$0.4d\n',...
   % 'Side stream flowrate,         S[kmol/min]      = %2$0.4d\n',...
    %'boilup flowrate,              V[kmol/min]      = %3$0.4d\n',...
   % 'reflux flow rate,             L[kmol/min]      = %4$0.4d\n',...
    %'liquid split ratio,           Rl               = %5$0.4d\n',...
   % 'vapor split ratio,            RV               = %6$0.4d\n',...
    %'impurity B in D,              xD,B[molB/mol]   = %7$0.4d\n',...
    %'impurity A in S,              xS,A[molA/mol]   = %8$0.4d\n',...
    %'impurity C in S,              xS,C[molC/mol]   = %9$0.4d\n',...
    %'impurity B in R,              xR,B[molB/mol]   = %10$0.4d\n',...
    %'composition A in feed,        zF1[molA/mol]    = %11$0.4d\n',...
    %'composition B in feed,        zF2,[molB/mol]   = %12$0.4d\n'...
   % ),u_nom(10),u_nom(5),u_nom(11),u_nom(2), u_nom(1),u_nom(3),u_nom(4),...
    %xf(im_xt(2)), xf(im_xS1(1)), 1-xf(im_xS1(1))-xf(im_xS1(2)),xf(im_xr(2)),...
    %u_nom(7),u_nom(8));
    
    petlyuk_cost
    %Calculate Juu
    

    %%
%Juu
%Juu = [1861 -1135 -17.9;-1135 1160 12.83;...
  % -17.9 12.83 1.84];                           when D and B loop was open
  %Juu = [1873 -1126 -17.72;-1126 1161 12.50;...      % value from Central diffrence close
 %-17.72 12.50 1.9609];
 Juu = [1858 -1131 -22;-1131 1145 8.12;...      % value from forward diffrence close
 -22 8.12 1.80];
 %Juu = [1857 -1131 -22;-1131 1147 8.12;...      % value from backward diffrence close
 %-22 8.12 1.81];
%Script for checking positive definite of a matrix
 %[m,n]=size(Juu); 
    %if m~=n,
     %  error('Juu is not Symmetric');
    %end
    
    %Test for positive definiteness
    %x=1; %Flag to check for positiveness
    %for i=1:m
       %subJuu=Juu(1:i,1:i); %Extract the upper left kxk submatrix
       %if(det(subJuu)<=0); %Check wheather the determinent of the kxk submatrix is +ve
          %  x=0;
           % break;
       % end
    %end
    
    %if x
       % display('Juu Matrix is Positive definite');
    %else
      %  display('Juu Matrix is NOT positive definite');
    %end      
