clc
clear all
close all

global h
h=actxserver('UnisimDesign.Application');

load x0_precooling
load dist_precooling

x0 = x0_precool ;

lb = [12/3600 ; 60.5/3600 ; 1100 ; 420 ; 210 ; 105 ] ;
ub = [13/3600 ; 64/3600 ; 1500 ; 480 ; 250 ; 125 ] ;

parameters.scale_x = eye(length(lb));
parameters.constraint = 'temperature' ;
% Vary flow of propane to NG precoolers and plot superheating of propane
% out of V102, both expressed in temperature and in enthalpy

% c3flow = linspace(lb(1),ub(1),30) ;
% sup_temp = zeros(size(c3flow));
% sup_enth = zeros(size(c3flow));
% 
% parameters.constraint = 'enthalpy' ;
% for i=1:length(c3flow)
%     x = x0 ;
%     x(1) = c3flow(i) ;
%     [C Ceq] = Precooling_constraints(x,d0,parameters) ;
%     sup_enth(i) = -C(1)+150 ;
% end
% figure(1),plot(c3flow*3600,sup_enth)
% 
% parameters.constraint = 'temperature' ;
% for i=1:length(c3flow)
%     x = x0 ;
%     x(1) = c3flow(i) ;
%     [C Ceq] = Precooling_constraints(x,d0,parameters) ;
%     sup_temp(i) = -C(1)+2 ;
% end
% figure(2),plot(c3flow*3600,sup_temp)

%Plot objective function against all variables!

NoOfPoints = 50 ;
X = zeros(6,NoOfPoints) ;
for i=1:length(lb)
    X(i,:) = linspace(lb(i),ub(i),NoOfPoints);
end

J = zeros(size(X)) ;
disp('Calculating objective function')
for i=1:length(lb)
    for j=1:NoOfPoints
    x = x0 ;
    x(i) = X(i,j) ;
    J(i,j) = Precooling_obj(x,d0,parameters) ;
    end
end 
disp('Done calculating objective function')
X(1:2,:) = X(1:2,:)*3600 ; % Change units for the two flows!

for i=1:6
    figure(i),plot(X(i,:),J(i,:))
end

%Repeat, for constraints. We already have the "X" matrix ready.

C1 = zeros(size(X)) ;
C2 = zeros(size(X)) ;
C3 = zeros(size(X)) ;
C4 = zeros(size(X)) ;
C5 = zeros(size(X)) ;
disp('Calculating constraint values')
for i=1:length(lb)
    for j=1:NoOfPoints
    x = x0 ;
    x(i) = X(i,j) ;
    C = Precooling_constraints(x,d0,parameters) ;
    C1(i,j) = C(1) ;
    C2(i,j) = C(2) ;
    C3(i,j) = C(3) ;
    C4(i,j) = C(4) ;
    C5(i,j) = C(5) ;
    end
end 
disp('Done calculating constraint values')
whattoplot = input('Which constraint should we plot? Type a number from 1 to 5: ') ;

switch whattoplot
    case(1)
        for i=1:6
            figure(i+6),plot(X(i,:),C1(i,:))
        end
    case(2)
        for i=1:6
            figure(i+6),plot(X(i,:),C2(i,:))
        end
    case(3)
        for i=1:6
            figure(i+6),plot(X(i,:),C3(i,:))
        end
    case(4)
        for i=1:6
            figure(i+6),plot(X(i,:),C4(i,:))
        end
    case(5)
        for i=1:6
            figure(i+6),plot(X(i,:),C5(i,:))
        end
end
    
disp('Resetting Unisim model')
Reset_Precooling(x0_precool,d0,h) ;