U = un(1:11);
nr = un(12:12+param.NC-1) ; 

disp(['Reactor temperature: ' num2str(U(1)) ' K']);
disp(['Reactor holdup: ' num2str(sum(nr)) ' mol']);
disp(['Reactor feed flowrate: ' num2str(U(2)) ' mol/s']);
disp(['Reactor feed molar fractions: ' num2str(U(3:5)')]);
disp(['Reactor exit flowrate: ' num2str(U(6)) ' mol/s']);
disp(['Reactor composition: ' num2str(nr'/sum(nr))]);
disp(['Column reflux rate (L): ' num2str(U(7))]) ;
disp(['Column boilup rate (V): ' num2str(U(8))]) ;
disp(['Fresh A feed rate (FA): ' num2str(U(9))]) ;
disp(['Purge ratio (P/D): ' num2str(U(10))]) ;


%Run "reactor.m once more to find reaction rates
Tr = U(1) ; %#ok<*NASGU> % Reactor temperature
RF = U(2)  ;% Flow to reactor
x_rf = U(3:5);  % Composition of reactor feed
RP = U(6) ; % Flow out of reactor
L = U(7)  ;% Column reflux rate
V = U(8)  ;% Column boilup rate
FA = U(9) ; % Fresh A feed rate
P = U(10)  ;% Purge ratio
RCY = U(11) ; % Recycle rate
reactor



%Calculate reactor yield/selectivity:
Y = convers ; 
S = 100*Y/conv_A ;

ncol = un((12+param.NC):end) ;
ucol = [L V] ;
dcol = [RP 1 xr'] ;
[dndt_col B D xB xD] = column(ncol,ucol,dcol,param) ;

disp(['Purge flow rate (P): ' num2str(P*D)]) ;
disp(['Recycle flow rate (RCY): ' num2str(U(11))]) ;
disp(['Generated by reaction: ' num2str(gen') ' mol/s'])
disp(['Column distillate flowrate: ' num2str(D) ' mol/s']);
disp(['Column distillate composition: ' num2str(xD')]) ;
disp(['Column bottoms flowrate: ' num2str(B) ' mol/s']);
disp(['Column bottoms composition: ' num2str(xB')]) ;

disp(['Conversion of A in reactor: ' num2str(conv_A) ' %'])


if ~isempty(param.prices)
    J = [FA V B (D-RCY)]*param.prices ;
    disp(['Value of objective function J: ' num2str(J) ])
end

NT = param.NT;
iu = 1:NT ;
xa = ncol(iu);
xb = ncol(NT+iu);
xc = 1 - xa - xb ;

figure(666),plot(iu,xa,'b',iu,xb,'r',iu,xc,'k')