% 10.34 - Fall 2006
% HW Set#2, Problem #4a
% Rob Ashcraft - Sept. 11, 2006

% This function takes in a vector of progress variables and reactor
% temperatures and calculates the residual values at each combination.
% Then it plots the residual values as a surface, as well as a contour plot
% showing the zero contour of both equations.  The intersection of the zero
% contours is the solution to the problem for the given conditions.

%%%%%%%%%%%%%%%%%%%%%%%%%
function WGS_contour(z_vec,Tr_vec,P,Tin,x_in)

global Vreactor res_time R Rj A Ea dHrxn dSrxn Cp_vec h

Cp_in = Cp_vec*x_in;


x_CO_in = x_in(1);
x_H2O_in = x_in(2);
x_CO2_in = x_in(3);
x_H2_in = x_in(4);

PCOin = P*x_CO_in;
PH2Oin = P*x_H2O_in;
PCO2in = P*x_CO2_in;
PH2in = P*x_H2_in;

h = 1.2; %1.2e-3;  %W/K

for i=1:length(z_vec)
    for j=1:length(Tr_vec)
        z(i,j) = z_vec(i);
        Treactor(i,j) = Tr_vec(j);
        
        dGrxn = dHrxn - Treactor(i,j)*dSrxn;
        
        PCO = PCOin - z(i,j)*R*Treactor(i,j)/Vreactor;
        PH2O = PH2Oin - z(i,j)*R*Treactor(i,j)/Vreactor;
        PCO2 = PCO2in + z(i,j)*R*Treactor(i,j)/Vreactor;
        PH2 = PH2in + z(i,j)*R*Treactor(i,j)/Vreactor;

        rate = A*exp(-Ea/Rj/Treactor(i,j))*(PCO*PH2O - PH2*PCO2*exp(dGrxn/Rj/Treactor(i,j)))/(1 + (PCO/0.2));

        resid_mass(i,j) = z(i,j) - rate*res_time;

        resid_energy(i,j) = -dHrxn*rate - (P*Vreactor/R/Treactor(i,j))/res_time*(Cp_in*(Treactor(i,j)-Tin)) - h*(Treactor(i,j)-300);
        
        
    end
    if(rem(i,100) == 0)
        disp(['Completed ',num2str(i),' steps out of ',num2str(length(z_vec)),' in making the contour plot']);
    end
end

z_max = x_in(1)*P*Vreactor/R/Tin;

z_mat = z;

Conv_mat = (z)/z_max*100;
T_mat = Treactor;

figure;
[C,h] = contour(Conv_mat,T_mat,resid_mass,[0 0],'-.k');
text_handle = clabel(C,h);
set(text_handle,'BackgroundColor',[1 1 .6],...
    'Edgecolor',[.7 .7 .7]);
hold 'on'
[C2,h2] = contour(Conv_mat,T_mat,resid_energy,[0 0],'--b');
text_handle = clabel(C2,h2);
set(text_handle,'BackgroundColor',[1 1 1],...
    'Edgecolor',[.7 .7 .7]);
legend('Mass Balance Residual','Energy Balance Residual');
xlabel('CO conversion (%)'); ylabel('Reactor Temperature (K)');
title(['Zero-Residual Contours of the Reaction and Energy Balances for T_I_n = ',num2str(Tin),' K']);


figure;
Cm=(zeros(length(z_vec),length(z_vec),3));
mesh(Conv_mat,T_mat,resid_energy,Cm);
xlabel('CO conversion (%)'); ylabel('Reactor Temperature (K)');
title(['Residual Values for the Energy Balances for T_I_n = ',num2str(Tin),' K']);
axis([0 1 0 1 -1e7 1e7]); axis 'auto xy'; 

figure;
surf(Conv_mat,T_mat,resid_energy,'EdgeColor','none');
xlabel('CO conversion (%)'); ylabel('Reactor Temperature (K)');
title(['Residual Values for the Energy Balances for T_I_n = ',num2str(Tin),' K']);
colormap('hsv');
axis([0 1 0 1 -1e7 1e7]); axis 'auto xy';  caxis([-1e7 1e7]); colorbar;

figure;
surf(Conv_mat,T_mat,resid_mass,'EdgeColor','none');
xlabel('CO conversion (%)'); ylabel('Reactor Temperature (K)');
title(['Residual Values for the Mass Balances for T_I_n = ',num2str(Tin),' K']);
colormap('hsv');
axis([0 1 0 1 -1e2 1e2]); axis 'auto xy';  caxis([-1e0 1e0]); colorbar;