%************************************************************** 
%        VARYING HEAT LOAD (Q_HP) TO THE HP COLUMN
% 
% - calls broyden1.m to find solution in one initial point and
%   to estimate the Jacobi matrix
% - then calles secant_cont.m 
% - inputs: Results_steady_state
% - plots various variables as a function of the heat load (Q_HP)
%**************************************************************
NT = 16;        % Number of stages in HP column
Nt = 26;        % Number of stages in LP column

load Results_steady_state;

u = [D_HP;B_HP;LT_HP;D_LP;B_LP;VD_LP;LT_LP;VB_HP];
ux0= [Xss;u];       % Initial values for broyden routine
B = [];             % Jacobi matrix

% Run the secant_cont program
%ux0 = ux;              % Initial solution vector
y = ux0(134);           % 'Free variable value' in last solution
yold = ux0(134)-ux0(134)*0.001;    % 'Free variable value' in solution *0.001
                                   %  before last solution
da = 1;                % Steplength for each iteration
mc = 4000;             % Number of iteration steps 
mi = 100;              % Maximum number of iterations
varargin = [];

[ux,B,step1,tt] = secant_cont('multieffect_heatload',ux0,y,yold,B,da,mc,mi,varargin);

clear B step1 tt
ux=ux([1 NT 2*NT+1 3*NT 3*NT+1 3*NT+Nt 3*NT+2*Nt+1 127 128 130 131 134],:);
% Results
xB_HP   = ux(1,:);    % Btm mole fn in HP column
xD_HP   = ux(2,:);    % Top mole fn in HP column
TB_HP   = ux(3,:);
TC_HP   = ux(4,:);    % Temperature in top HP col (K)
xB_LP   = ux(5,:);    % Btm mole fn in LP column
xD_LP   = ux(6,:);    % Top mole fn in LP column
TB_LP   = ux(7,:);    % Temperature in btm LP col (K)
D_HP    = ux(8,:);    % Distillate flow in HP column
B_HP    = ux(9,:);    % Bottom flow in HP column
D_LP    = ux(10,:);    % Distillate flow in LP column
B_LP    = ux(11,:);   % Bottom flow in LP column
VB_HP   = ux(12,:);   % Hvap in HP reboiler (kJ/kmol)

% Calculate pressure on top stage in HP column
R = 8.314;                     % Universal gas constant (kJ/kmol*K)
TBL = 385.85;TBH = 426.45;     % Boiling point pure components (K)
HvapL = 31630; HvapH = 37870;  % Hvap for pure components (kJ/kmol)
p0L = 5.2e5; p0H = 5.2e5;      % Vap. pressure for pure liq. comp.(Pa)
pL = p0L*exp(-HvapL/R*(1./TC_HP - 1/TBL));  % Vap. pres. light comp.
pH = p0H*exp(-HvapH/R*(1./TC_HP - 1/TBH));  % Vap. pres. heavy comp.
Ptop_HP =  xD_HP.*pL + (1-xD_HP).*pH;       % Stage pres.(Raoult's law)

% Heat duty in HP column
Hvap_btm_HP = 35440;         % From HYSYS at 5.2bar & 0.53 MeOH 
Qb_HP   = Hvap_btm_HP*VB_HP/(3600*1000);       % kJ/h --> MJ/s

% Objective function (mol/s)
wr = 0.6488;                   % mole/MJ
J  = (D_HP+D_LP).*(1000/3600) - Qb_HP*0.6488;    % Objective function (mol/s)


Results = [D_HP' B_HP' xD_HP' xB_HP' D_LP' B_LP' xD_LP' xB_LP' TB_LP' TC_HP' VB_HP' Qb_HP'  Ptop_HP' J'];

% Save results
% save Results_vary_heatload Results 

% Plot Results
figure(1); 
subplot(3,1,1); plotyy(Results(:,12),(1-Results(:,8)),Results(:,12),Results(:,9));
title('Temperature and Water Composition in LP Column')
ylabel('Water conc. btm. LP col.')
subplot(3,1,2); plot(Results(:,12),Results(:,13)/1e5);
title('Pressure (at top) in HP Column')
ylabel('Pressure (bar)')
subplot(3,1,3); plotyy(Results(:,12),(1-Results(:,4)),Results(:,12),Results(:,2));
title('Bottom Flow and Water Composition in HP Column')
xlabel('Heat load (MW)'); ylabel('Water conc. btm HP column')
hold

figure(2); 
subplot(2,1,1); plot(1-Results(:,4),Results(:,14));
title('Objective function vs Btm. Water Composition in HP Column')
xlabel('Btm. water composition'); ylabel('Objective')
subplot(2,1,2); plot(Results(:,2),Results(:,14));
title('Objective function vs Btm. Flow in HP Column')
xlabel('Btm. flow, B_{HP} (kmol/h)'); ylabel('Objective')
hold

figure(3);
subplot(2,1,1); plot(Results(:,12),(Results(:,1)+Results(:,5))*(1000/3600));
title('Production Rate as a fn of Heat Load')
xlabel('Heat load (MW)'); ylabel('Produced methanol (mol/h)')
subplot(2,1,2); plot(Results(:,12),Results(:,14));
title('Objective function as a fn of Heat Load')
xlabel('Heat load (MW)'); ylabel('Objective function')
hold



