%% Published by Bjorn Tore Mathisen
% Model of the Bulk Cooler in NPK 4 as part of a specialization 
% subject H2013.
% The following model approximates the Bulk Cooler as a CSTR:

clc
clf
clear all
close all

% Process constants:
constant

% Global parameters
parameters

%%
% Disturbances
% time_n_R_BC = 20;        % Steptime [min]
% step_n_R_BC = 0.01;         % Stepfac from init value
% time_n_BC_S = 20.5;        % Steptime [min]
% step_n_BC_S = 0.01;         % Stepfac from init value
time_T_R      = 20;
step_T_R      = 0.0;

%Starting simulink model
sim('CSTRseries.mdl');


disp('Total energy balance:')
text1 = ['Qwater: ',num2str(i.n_WR_W*4183*(x(end,1+i.n+3)-c.T_cin)),' J/s'];
disp(text1)
text2 = ['QNPK: ',num2str(i.n_BC_S*p.a_NPK1*(c.T_R-x(end,10+3))),' J/s'];
disp(text2)

figure(1)
subplot(211)
hold on
stages1 = [1:5];
axis([0 6 12 45])
title('Top section: Counter current Hx')
plot(stages1,x(end,stages1+2),'-xr')
plot([0 1],[c.T_R x(end,3)],'-xr')
plot(stages1,x(end,stages1+i.n+3),'-sb')
plot([5 6],[ x(end,5+i.n+3) x(end,10+i.n+3)],'-sb')
plot([0.5 0.5], [13 44],'--k')
plot([5.5 5.5], [13 44],'--k')
str(1) = {'T_{R}'};
text(0.1,38,str(1))
str(2) = {'Inlet'};
text(0.1,(45-12)/2+12,str(2))
str(6) = {'Inlet'};
text(5.6,(45-12)/2+12,str(6))
str(3) = {'Hx'};
text(0.6,(45-12)/2+12,str(3))
str(3) = {'Hx'};
text(5.1,(45-12)/2+12,str(3))
str(4) = {'T_{10}'};
text(5.6,16,str(4))
ylabel( 'T [C]')

figure(1)
subplot(212)
stages2 = [6:10];
hold on
axis([5 11 12 45])
title('Bottom section: Co-current Hx')
plot(stages2,x(end,stages2+3),'-xr')
plot( [5 6],[x(end,5+2) x(end,6+3)],'-xr')
plot(stages2,x(end,stages2+i.n+3),'-sb')
plot([5 6],[c.T_cin x(end,6+i.n+3)],'-sb')
plot([5.5 5.5], [13 44],'--k')
str(1) = {'T_{cin}'};
text(5.1,18,str(1))
str(2) = {'Inlet'};
text(5.1,(45-12)/2+12,str(2))
str(3) = {'Hx'};
text(5.6,(45-12)/2+12,str(3))
str(4) = {'T_{5}'};
text(5.1,38,str(4))
ylabel( 'T [C]')
xlabel('Stage')

% figure(2)
% subplot(211)
% stages2 = [1:12];
% hold on
% axis([1 12 15 45])
% plot(stages2,x(end,stages2+1),'sr')
% ylabel( 'T [C]')
% xlabel('Stage')
% 
% figure(2)
% subplot(212)
% stages2 = [6:10];
% stages1 = [1:5]
% hold on
% axis([1 10 15 45])
% temp = [x(end,stages2+i.n+3) x(end,stages1+i.n+3)];
% plot([stages1 10 9 8 7 6 ],temp,'sb')
% ylabel( 'T [C]')
% xlabel('Stage')